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FOREWORD 


This interim report, "Controller Design Technology for the Space Shuttle 
Vehicle; summarizes the study performed during the period 10 June 1970 
to 1 June 1971 for the National Aeronautics and Space Administration, 
George C. Marshall Space Flight Center, under Contract NAS8-25708. 

Dr. S. W. Winder of the Dynamics and Control Division of the Aero- 
Astro'dynamics Laboratory was the technical monitor . The study was per- 
formed in the Systems and Research Division of Honeywell Inc. Dr. G. B. 
Skelton served as program manager. Dr. A. J. VanDierendonck was the 
initial principal investigator. He completed the study of iteration methods 
for the time -varying problem. Dr. C. A. Harvey, succeeded Dr. 
VanDierendonck as principal investigator, completed the work on controller 
simplification and conducted the study concerned with parameter variations. 
Dr. G. Stein became a coprincipal investigator and conducted the sensor 
choice study. Dr. Y. 'S. Lee contributed to the investigation of sensitivity 
to parameter variations. Mr. M. D. Ward assisted with computer pro- 
gramming and numerical analysis. 
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SECTION I 

INTRODUCTION AND SUMMARY 


The goal of the research program summarized in this report was to develop 
a technology for control system design. The technology was aimed at pro- 
viding control of flexure and rigid-body degrees-of -freedom of the Space 
Shuttle V ehi'cle within constraints of practicality. The constraints apply to 
the controller configuration and to the design methods as well. 

At the time the program was initiated, it was anticipated that the Space 
Shuttle Vehicle would exhibit flexure control problems caused by large sizes 
of booster and orbiter sections, coupled modes of the vehicles in mated 
ascent, sensitivity to gust and maneuver excitation caused by large aero- 
dynamic surfaces, and possible orbiter structural requirements imposed 
by reentry heating.. The large aerodynamic surfaces could also cause a 
significant load relief problem at maximum dynamic pressure and larger 
drift dispersions than were present with Saturn Vehicles. Honeywell had 
achieved considerable success in treating such problems for large launch 
vehicles and large flexible aircraft. A stochastic constrained -response 
theory was developed, and its applicability to the control of a rigid booster 
was demonstrated in 1965 and 1966 (Ref. 1). The B-52 LAMS (Load Alle- 
viation and Mode Stabilization) system was designed and flight-tested in 
1966 and 1967 (Ref. 2). In 1967 and 1968, applicability of the method to a 
flexible launch vehicle was demonstrated (Ref. 3). In spite of the success 
achieved, the technology developed by 1968 was inadequate with respect to 
three control design problems: (1) controller simplification, (2) sensitivity 
to model inaccuracy, and (3) sensor complement choice. 
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Thus the specific goals of this study were to improve the design technology 
by providing: 

• Practical controller simplification algorithms 

• A mathematical method for implicitly including parameter 
variation constraints within quadratic optimization formulations 

• A rigorous mathematical basis for understanding best sensor 
choice and location 

The first objective is motivated by the difficulty of past efforts to simplify con- 
trollers. Simplification of the optimal controller in the LAMS program required 
a nine -man /month simulation program, and the optimal controller for the 
flexible launch vehicle was reduced by similar expensive trial and error 
methods. At that time, necessary conditions for optimizing constrained 
controller configurations were known (Ref. 4). These conditions provided a 
two-point boundary -value problem and computational algorithms for its 
solution. However, it was too expensive to use these algorithms for such 
problems as the flexible launch vehicle control problem, especially if 
several measurement complements or sets of parameter values were to be 
considered. Another drawback to these algorithms was the possibility of 
convergence to an arbitrary local minimum. Improved computational pro- 
cedures for the solution of optimal control problems and simplification of 
optimal controllers had been developed at Honeywell by 1970 which held 
promise for treating the control problems of highly flexible vehicles. During 
the present study, these procedures were used as a basis for developing 
practical controller simplification algorithms. For the Space Shuttle Vehicle 
these algorithms may be used for controller design for mission phases such 
as the highly flexible orbiter during re-entry and cruise or for control during 
the mated ascent phase. 
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Mathematical models of two test vehicles were used for assessing the 
capability of the simplification procedures and improving them where 
necessary. A rigid representation of the Saturn V with a Voyager payload 
was used for a time-varying ascent study, and a model of the longitudinal 
axis of the C-5A with six flexure modes at a single flight condition was used 
for study of flexible -vehicle control. These choices were made to reduce 
the development cost, but the results of the two studies could be combined for 
actual application to Space Shuttle Vehicle controller design. 

Controller simplification in both cases started with the solution of the 
stochastic constrained-response optimization problem with complete state 
measurement capability. For the launch vehicle, this solution was a set of 
time -varying gains and a time -varying deterministic input which defined the 
controller. This controller minimized a quadratic performance functional 
and, via quadratic equivalence, minimized as well a nonquadratic performance 
functional called the cost functional. This cost functional represents an upper 
bound on the likelihood of mission failure. A gradient iteration technique for 
controller simplification along with methods for choosing initial conditions 
for the iteration technique were developed and tested. These techniques 
proved to be quite capable for the example treated. The original controller 
utilizing ten time -varying feedbacks was simplified to a controller with five 
time -varying feedbacks. 

Attempts at controller simplification for the flexible C-5A vehicle with an •, 
Implicit Function Method were unsuccessful. The cause of the failure was 
attributed either to lack of damping in the algorithm used or to severe 
sensitivity of performance to gain changes. Damping was added to the 
algorithm. This yielded some improvement, but not enough to provide a 
solution within assumed practicality constraints. Extrapolation based on the 
initial step of the Implicit Function Method was then used with gradient 
correction to achieve a satisfactory solution. A. third method was also used 
successfully. This method, called the Incremental Gradient Method, 
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incorporates desirable features of the Implicit Function Method and usual 
gradient techniques. 

The description of both phases of controller simplification is presented in 
Section II and Appendices A through C. 

Another area in which the design technology needed upgrading to be applicable 
to the Space Shuttle Vehicle was sensitivity to model inaccuracy. In every 
iteration of controller design there is a degree of uncertainty in the data. 

Also, certain dynamics such as high order flexure and fuel- sloshing modes 
are generally ignored to make the design problem tractable. Thus, practical 
controller design must recognize and be tolerant of model inaccuracies. In 
this study such inaccuracies were assumed to take the form of unknown param- 
eters. Two directions were pursued to yield formulations of the quadratic 
optimal control problems which would implicitly include parameter uncertainty 
constraints. One direction was to determine properties of the optimal per- 
formance surface over a segment of parameter space which could be used to 
-derive optimal insensitive controllers. An optimal insensitive controller is 
'one that is optimal for a given value of the parameters defining the system 
and minimizes the maximum of the performance index over all admissible values 
of the parameters. A necessary condition for an optimal insensitive con- 
troller was derived. This condition was shown to be locally sufficient. The 
nature of these conditions led to the development of an algorithm for computing 
approximately optimal insensitive controllers. The utility of this computational 
approach was tested on the C-5AT example, and significant reduction in sen- 
sitivity was achieved. 

The other direction of the parameter variation study was an investigation of 
the effectiveness of compensators in reducing sensitivity. This study led to 
choosing control parameters to match the performance of dominant dynamics 
of the compensated system to the performance of the optimal uncompensated 
system over the range of system parameters. This approach is intuitively 
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appealing. But for high -dimensional systems the approach is of questionable 
value. 

Analyses performed in these methods of including parameter variation con- 
straints are described in Section III and Appendix D. Results derived for 
approximately optimal insensitive controllers for the C-5A example are 
presented in Appendix E. 

A third area of technology improvement concerns the choice of types and 
locations of sensors to generate feedback signals for practical control systems. 
This so-called "sensor -choice problem" arises because it is economically 
prohibitive in most applications to measure all system states, particularly 
the rates and displacements of flexure modes. The problem is unsolved 
because it is theoretically and computationally difficult to determine basic 
performance capabilities of a set of sensors. To do so requires the solution 
of two coupled optimization problems — (1) optimal location of the instruments 
and (2) optimal design of a practical controller to utilize the instruments. 

Both are complex problems when a quadratic performance functional is used 
as the measure of quality. Available methods are based on controller 
optimization routines such as those discussed in Section II nested within 
iterative search procedures for best sensor locations (Ref. 5, 11). 

In this study, the sensor choice problem was approached from the viewpoint 
of finding alternate quality measures which would be more convenient com- 
putationally yet still provide meaningful indications of performance capability. 
Two measures were proposed, both based on the pole -placement capability 
of the sensor complement. The first measure is the maximum number of 
closed-loop poles, p , which can be placed arbitrarily, and the second is a 
measure of deviation of the remaining nonarbitrary poles from specified 
desirable locations. General equations were derived for these quantities, 
and a computational pole placement algorithm was developed for their solution. 
Computational feasibility of the measures was established. However, their 
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practical utility as quality measures, particularly their correlation with 
quadratic cost, remain to b,e verified. 

Analyses performed on the sensor choice problem .are described in Section IV 
and Appendices F, G .and H. 
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SECTION II 

CONTROLLER SIMPLIFICATION 


The aim of this phase of the study was to develop economically feasible 
methods for computing simplified controllers for the Space Shuttle Vehicle. 
The simplified, controllers were to be derived from optimal controllers. 
Simplification was defined to be the reduction in the output measurements 
required and, in the case of time-varying gains, reduction in the complexity 
of the time-variations to specified parametric representations. Of course,’ 
the simplified controllers were to maintain as much as possible the desired 
performance of the optimal controllers. 


BACKGROUND 

The mission of the Space Shuttle Vehicle is such that at times the control 

\ 

problems are similar to those of a large flexible launch vehicle and at ^ 
other times to those of large flexible aircraft. The stochastic constrained- 
response formulations of these two types of control problems had been 
derived in References 1 and 2. Resulting optimal controllers had been shown 
to provide very desirable performance. The measures of performance and 
the mathematical models for large flexible launch vehicles and large flexible 
aircraft have several common features. The major distinction from a 
mathematical viewpoint is that the control problem for the launch vehicle is 
a finite- time problem with significantly time-varying dynamics, while the 
aircraft can be considered to fly at a single flight condition for a long enough 
period of time that the dynamics maybe represented by a constant-coefficient 
model, and a steady-state .performance functional is of interest. Common 
features include linear dynamics, many degrees of freedom, stochastic 
disturbance models, and physically meaningful performance functionals 
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which are generally not expressable as integrals of quadratic functions of 
(linear) responses. These common features suggest that, in theory, 
simplification methods which are successful for one problem will also be 
successful for the other. In practice, the computational costs for the two 
problems differ significantly. For example, with a quadratic performance 
functional, the cost of computing an optimal launch vehicle controller for- a 
23rd-order Saturn model is approximately 50 times the cost of computing 
an optimal controller for a 2 3rd- order model of a flexible aircraft at a 
single flight condition. These facts motivated the separation of the con- 
troller simplification study into two parts, one part dealing with the time- 
varying aspect of the problem and the other dealing with the high dimensionality 
of the mathematical models associated with flexure. 

In addition to the difference in computation costs between the time-varying 

and constant- coefficient problems, there is a difference in computer storage 

requirements. For problems with the same number of degrees of freedom, 

the time-varying problem requires much, more computer storage than does 

the constant-coefficient problem. Computation time considerations led to 

the choice of a rigid model of a launch vehicle as a test vehicle for the time- 

varying aspects of this study. This permitted savings in computer costs 

* 

because of the reduction in order of the system from the model with flexure 
and also because larger sampling intervals could be used. The flexure 
degrees of freedom were retained in the constant-coefficient model. 

Computation time and storage requirements also influenced the type of 
techniques chosen for simplification in the two problems. The high cost of 
computation for the time-varying problem limited the acceptable m'ethods 
to first-order, that is, methods which utilize at most first-order derivatives. 
Large storage requirements impose the same limit for the. constant- 
coefficient problem if the number of nonzero parameters (gains) permitted 
in the simplified configuration is large. If the number of such parameters is 
not large, then second-order methods are candidates for consideration. 
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SIMPLIFICATION PROCEDURES 


The starting point for the procedures is a solution to the following 
optimization problem. Dynamics of the system are represented by the 
linear differential equation 

x = Fx + G 1 u-+G 2 r| (1) 

where x is the state vector, u is the control vector and r) is a white-noise 
vector of disturbances with known means and covariance. Responses to be 
controlled are linear combinations of the states and controls given by 

r = Hx + Du (2) 

The performance index is a functional of the form 

\ 

J = f 1 [R(T>,S(T)]+ J' T f 2 [R(t) 3 S(t)]dt (3) 

o 


where 


R(t) = r(t)[r(t)l T , r(t) = E{r(t)} 

(4) 

S(t) = E{ [r(t)-r(t)l [r(t)-r(t)] T | 

(5) 


The superscript T, denotes the transpose. The optimization problem is to 
choose u to minimize J. For the ascent problem, the time T in (3) is finite 
and F, G^, G 2 , H, D, the function f^, and the mean and covariance of the 
noise maybe explicit functions of time. For the flexible- aircraft problem, 
the coefficient matrices, the function f 2 , and the noise characteristics are 
constant, f^ is zero and T in the integral in (3) is infinite. 
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The "quadratic problem" occurs when and f 2 in (3) are linear functions 
of K and S, since these matrices are quadratic in the responses. The 
solution of the quadratic problem is 

u = K*x + f (6) 

where K* and f satisfy certain determining equations. For simplicity of 
discussion, let us consider the constant- coefficient problem with 
f = Tr[QSl where Tr denotes the trace. Also, for simplicity assume 
E |T) | = 0 which implies f = 0. This last assumption does not represent a 
restriction on the simplification problem. Inclusion of a deterministic input 
in the controller is not a real complication. The determining equations for 
K are 

0 = D T Q(H+ DK*) + G 1 T P ■ '(7) 


and 

0 - (F + GK*) T P + P(F + GK*) + (H + DK*) T Q(H + DK*) (8) 

Analytically, controller simplification consists of replacing the matrix K* 
with a matrix K of a constrained form which has less independent parameters 
than the number of elements in K*. The resulting controller, of course, should 
maintain as low a value of J as possible. Thus at this point, J may be con- 
sidered as a function of K to be minimized subject to the constraining con- 
figuration imposed on K. The procedures used in this study for performing 
this task may be classified as gradient iteration and parametric techniques. 


Gradient Iteration Technique 

The gradient iteration technique consists of applying gradient methods to the 
minimization with respect to K of the Hamiltonian 
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H(X, P, K) = Tr/Q(H + DK)X(H + DK) T + P[(F + G 1 K)X+X{F + G 1 K) T +N](t 

(9) 

rp rn 

where X denotes the state covariance matrix and N = EjGgfM G 2 j.. 

The necessary conditions for optimality of K derived by Axsater (Ref. 4) are 
d H . ?H 3H 

X = — , p = , — = 0 (10) 

sp dX 3K 

with K constrained to be of the desired configuration. If complete measure- 
ment of the state is permitted, the first of equations (10) is uncoupled from 
the last two. With complete measurement of the state not permitted, the 
three equations are coupled. Thus at each iteration, the first two equations 
are integrated to find X and P. Then the gradient 9H/9K is evaluated. When 
this gradient is nonzero, a step is taken to reduce the value of the Hamiltonian. 

Existence of several local minima for the constrained minimization problem 
is indeed possible. But if an initial value for K is known which is sufficiently 
close to the global minimum, the gradient methods have the desired property 
of convergence to the global minimum with sufficiently small steps. Thus, 
the successful application of this technique depends on having a good initial 
value for K, proper choice of step size, and, of course, the capability to 
compute gradients economically. 

Two methods were used to obtain initial values of K for the iteration studies. 
The assumed configuration for the controller was u = Kz where 

z = Mx (11) 

tThe second term in the trace is identically zero for the s'teady-state problem. 
It is included here to indicate the presence of its analog in the time-varying 
problem. 
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represented a vector of measurements. The first method was to choose K 
to be the projection of K* on M 

K = K*M T (MM T )' 1 (12) 

When the components of z are components, of x, the choice in (12) merely 
deletes the gains in K* which multiply states that are not measured. The 
second method was based on approximating a Kalman Estimator used in 
conjunction with K* by the gain K. 

Essentially, this method involves finding a good. match between the open-loop 
responses of Figures 1(a) and 1(b), for a very small measurement noise S. 
Figure 1(a) yields a response in the frequency domain of 

U(jw) = Y(j<u)Z(ju) (13) 

where Y(s) is the Laplace transform of the system shown. Specifically 

Y(s) = K(sl - F - G-jK* + LM) -1 L (.14) 

where L is' the Kalman gain matrix. 

Figure 1(b) yields a response of 

U(jw) = KZ(j w) (15) 

where K is not a, function of frequency. 

Two methods' to match these responses were tried. The- first was to set the 
static gains equal to the d-c gains of the dynamic, system.- That is 

K =' Y(0) (16) 
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CONTROL INPUT TO KALMAN ESTIMATOR 



? (MEASUREMENT NOISE) 


(a) DYNAMIC 



(b) STATIC 


Figure 1. Dynamic and Static Gain Systems 

The second was to compute a weighted average of Y(j^) over all frequencies. 
That is 

CO 

Re jY(ju>)|Z(ju>) |d«J 

K = (17) 

do 

J |Z(jw) |d«J 

o 

where Re denotes "the real part of. " These methods were actually applied 
to the time-varying example discussed below. 

The choice of step size is a problem common to all gradient iteration methods. 
Furthermore, it is well-known that the path of steepest descent is dependent 
on the coordinate system chosen. For simplicity of exposition, consider K to 
be a vector with components k.. The gradient iteration maybe expressed as 
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K(p+1) = K(p)-eV K H 


e >o 


(18) 


| K=K(p) 

where p denotes the stage of the iteration, and e denotes the step size. A 
sketch of possible behavior of such an iteration method is shown in Figure 2 
for a two-dimensional K. In this sketch, points A, B, C, D and E indicate 
possible values of K(p+1) corresponding to different choices of e. If the value 
of e chosen is less than or equal to the value corresponding to C, the gradient 
at K(p+1) would be directed in the general direction of the point where H = 0. 

A common manner of choosing e is to perform the minimization of HCK(p+l)] 
with respect to -the parameter e. The result for K(p+1) on the segment from 
K(p) .to C is depicted as point B. F.or problems such as the time-varying 
control problem, even this one dimensional minimization can be costly so 
that an approximation is generally accepted. It is also clear from Figure 2 
that admissible directions of the step from K(p) to yield decreases in Hare 
all directions with positive projections on the vector from K(p) to B. In 
the example treated in this study .such a modification was introduced. 

Equation (18) was 'replaced by 



-Figure 2. Gradient Step Possibilities 
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K(p+1) = K(p) - e(diaglk.(p)| ) V k H[ |v k h| l" 1 , e >0 (19) 

*- -*K=K(p) 

where diag j k^(p) j denotes the diagonal matrix with fk^(p) | as its i^ 1 diagonal 
element. In this form e represents the maximum percentage change in any 
element of K(p). This iteration equation proved to be successful in the 
example. 

The equations for computation of the gradient in the example are described in 
the discussion of the example. The required computation is the solution of 
a covariance and a co-state equation. For the steady-state problem, the 
gradient of H can be computed as follows. The steady-state covariance 
matrix, X, satisfies 

(F + G 1 K)X + X(F + G 1 K) T + N = 0 (20) 

Thus, H = Tr(H +DK) T Q(H + DK)X = J. The adjoint matrix A corresponding 
to X satisfies 

(F + G 1 K) T A + A(F + G 1 K) + (H + DK) T Q(H + DK) = 0 (21) 

Then 

9 J rp ; ^ rn ° X 

= Tr [2 (H+ DK) Q DE J X+ (H+ DK) Q (H+ DK) 3 (22) 

dK.. 3K, . 

D D 

. . A dK 3X 

where E 1 ^ = and is defined by 

dK. . SK. . 

1 D i] 
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( 23 ) 


(24) 

Thus, to evaluate all of the first partials of J, we need only evaluate the two 
covariance- type -equations (20) and (21) and the algebraic equations (24). This 
convenient property was first discovered by T. L. Johnson (Ref. 5). 

Parametric Techniques 

In this technique developed by G. Stein (Ref. 6), a general .gain matrix K is 
partitioned into two orthogonal -components: 

K = (K 1 + K S )M + XK 2 (25) 

where (K 1 ■+ K 3 )M(K 2 ) T = 0. 

The first component (K 1 + K 3 )M satisfies the constraint and hence represents 

-2 

gains to be retained. The second component X-K consists of gains to be 
discarded. The scalar 'parameter, X, is introduced in equation (25) to permit 
discarding the second component while maintaining one constraint on the first 

component. For generality, the first component is divided into two orthogonal 

1 3 1 3 

components EX and.K where K represents gains to be optimized and K is 

fixed. The necessary condition for optimality of with respect to the con- 
strained problem is 


rp i.: ii t 

(F + G K) + (F + G-K) 1 -+-G-E’ 13 X + X (G. E 3 ) = 0. 

1 aK.. -3K.. 
i] il 

Using the adjoint property of X and A 

= Tr .[2 <(H + DK) T Q DE* 3 X+.A (G^ 13 X + X <G 1 'E 13 ) T ] 

3K ii 

= 2 Tr j [(H + .DK) T QD + AG 1 3 E 13 x} 
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(26) 


^[(K 1 + K 3 )M + XK 2 ] 


3K 


1 


A = 0 


= 0 


Solutions to this equation may be obtained as noted by D. K. Scharmack 
(Ref. 7) by starting with a K in (25) with X = 1 that satisfies 


aj(K) 

= 0 


3K 


1 


(27) 


and choosing K X as a function of X such that equation (27) is maintained for 
0 ^ X=£ 1. The most appealing choice for X = 1 is the K that minimizes J 
globally; namely 

K* = (K 1 (l) +K 3 )M + K 2 (28) 


According to the Implicit Function Theorem (Ref. 8), K 1 (X) satisfying (27) is 
defined by 


dKV) 


dX 


a^UK 1 + k 3 )m + xk 2 ] 


i -l 


1 IT 
9K1K 


3 2 J[K X + K 3 )M + XK 2 ] 


(2 9) t 


Some computational difficulties were experienced in applying this technique 
using numerical integration of (29) for the C-5A example. It was noted that 
equation (29) is equivalent to 


1 

tEquations (26), (27) and (29) make sense in vector matrix notation if K x 

is written as a column vector, which is assumed. 
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d 

— (VJ) = 0 

dX 


( 30 ) 


where VJ is abbreviated notation for the partial of J with respect to the gains 
in K 1 . As a differential equation in X, equation (30) is neutrally stable, so 
a damping term of a(X) J was added to the right-hand side. The scalar a(^) 
was chosen to be positive since the integration is in the negative X direction. 
This modification gave some improvement, but not complete^success, A 
hybrid technique using this method of finding K (. 95) from K (1), then 
extrapolating to K X (0) and then using gradient correction was then used 
successfully. 

A variation on this theme was developed by A. J. Van Dierendonck (Ref. 9). 
This variation is called the incremental gradient algorithm. This algorithm 
begins with the same values of K and K 1 at A = 1. Then x is incremented to 
X + A X (with -1 < AX< 0). Several gradient steps are then taken to obtain 
K 1 (X + AX). Another increment in X is then taken with linear prediction 
used for K 1 . Again gradient corrections in K 1 are made and the process 
continued until X = 0. This method was also tried successfully on the C-5A 
example. 

A comparison of the three parametric methods is illustrated for a simple^ 
example in Figure 3. For ease of representation, it is assumed K is a two- 
vector and K 1 = k 1# K 2 = k 2 . In realistic examples, each of these axes may 
be multidimensional subspaces. The sketch is completely hypothetical. 

The Implicit Function Solution shown is assumed to be exact. The other two 
paths are intended to be only indicative of the actual situations. 
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< ^2 ^3 < ^4 < ^5 


Figure 3. Geometric Interpretation of Parametric Techniques 


19 


SIMPLIFICATION AND QUADRATIC EQUIVALENCE 

In the study of the time-varying example, it was discovered that the quadratic 
approximation to a nonquadratic performance functional used in quadratic 
equivalence is dependent on the constraints imposed. The following example 
demonstrates this phenomenon. 

Consider minimizing the function f(x, y) = x+y+3y for (x, y) belonging to 
certain sets, say K ± and Kg. We may think of these sets as representing 
sets of attainability for a dynamical system and suppose that and Kg 
correspond to complete and- incomplete measurements respectively. Then 
since measurements could be ignored in the complete measurement case, 

K 2 would be contained in K^ in fact, let us assume Kg is a proper subset 
of K^. To give this problem the interpretation of a stochastic optimal 
control problem, we may suppose that x and y represent covariances of two 
responses. This would imply x s 0 and v ^ 0. Now suppose that 

K^ = {(x,.y) : x^O and xy^l/3} 

Kg = {(x, y) : x^O and- xy&l} 

The minimum of f(x, y) for (x, y) in K^ occurs at x = 1, y = 1/3. This point 
also yields the minimum for the "quadratic" function f(x, y; q 1 , q 2 > = 
q x + q 2 y with q 1 = 1, q 2 = 3. Note that the vector (q^ q 2 > = (1, 3) is a 
normal to the boundary of K 1 at (x, y) = (1, 1/3). To interpret f(x, y; q^ <^) 
as a. "quadratic" cost function, we consider x and y as variances of 
responses V 1 and Vg so that f(x, y; q 1 , q 2 ) is quadratic in the responses. 

For (x, y) belonging to Kg the minimum of f(x, y) occurs at x = 2, y=l/2 
and normal at that point is (1, 4). Thus the equivalent quadratic for K^ is 
f(x, y; 1, 3) whereas for Kg it is f(x, y; 1, 4). Figure 4 shows the relevant 
geometry for this example. 
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This phenomenon may- be quite significant for problems' in which quadratic 
equivalence and simplification are considered.- The simplification algorithms 
for such problems should probably utilize" the gradients of the nonquad’ratic 
performance functional,. at" least -in- the final stages' and not rely on the 
quadratic approximation throughout the simplification procedure. This was' 
clearly the' cas'e id the study' of the following example. 


TIME^VAICYING EXAMPLE 


The vehicle- used'-in this example is the rigid-body representation of the- 
vehicle, described in Reference 3. The resulting model is 10th order and' 
includes, rigid body,'.. wind- filter, distributed wind loads', and gimbal” angle; 
For this model- of-' the- Vehicle, a sampling rate of 10 samples per second- 
was* adequate, rather than th'e 50 samples' per second required for the 
original-model which' included' three flexure 'modes. . 


Initial optimization' r-uris were' made’ as suming’-complete' measurement capa-' 
bility (i. e'.., every vehicle-'and'wirid state can be measured). Using-the con 1 
cept of-quadra v t ic - equivalence-, three- iterations yielded a controller whose-'" 
upper bound' of’ the^probability that mission failure (as- described in Reference' 3 
would-.occur was-- 1-.- 7*x 1-0 . The -upper bound on the- probability of mission 

failur.eds denotedby. J*, while denotes a'quadratic cost functional. ’ It'w'as- 
assumed-’ that 0 (pitch), 0' (pitch rate), *z (drift), z (drift rate), and |3- (gimbal 
angle)- could betmeasure d 'di-re'cfly but’ that the measurements were' noisy; ' A" 


Kalman JSstimdtordvas derived 10 'ex-tima'te u and x (win'd states-) arid xl, x„-, 

1 . 

and’ Xg. (load- distribution -states')-. • The- '-costs,- J* and ; 'J** were computed'-with 
the- ’optimum' gains and- -the'tECai-m'ah ’Estimator. Although the' quadratic-cost' 

increased-. slightly' over-' the perfect- sensing' case, the upper bound' on -the 

- , , , . , - — 4 " * - * # *' * ' 

probability of mis s ion f a'ilu r e lower e d : .'signif ic ah'tly from 0. 7 x 10 to '0; 3 x' 

10' 6 . 
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The desired form of simplified controller was chosen to be feedback of the 
measured quantities without using compensators. The simplest choice of 
initial gains for the gradient iteration method is the perfect sensing gains 
for states which can be measured and zero for the states not measured. 

For this example, this choice produced results comparable to the perfect- 
sensing gains. This was understandable since, for this rigid-body model, 
essentially all costs were due to the second bending moment at maximum 
dynamic pressure, and this can be controlled with pitch, pitch rate, and 
the gimbal angle. These results also explain why the Kalman Estimator 
did not appreciably change the results from those for perfect sensing. 

To create a problem to test the gradient iteration method, the terminal- 
drift constraint was reduced by a factor of three from 30, 000 to 10, 000 
meters. It was found that this new constraint was easily met. Hence, this 
did not create a meaningful problem. 

Thus, another approach was taken. Instead of tightening the constraints, 
the magnitude of the disturbance input was increased by doubling the 
standard deviation of the wind. After changing the quadratic weights, a 
new optimal controller was computed that met the original constraints. This 
controller resulted in an upper bound on the probability of mission failure 

_ r>o 

(J*) of 0. 4299 x 10 . The derivatives of this probability with respect to 

the individual responses indicated that quadratic equivalence was very 
nearly achieved. 

The value of J* for the controller with perfect-sensing gains for measured • 

— fi 

states and zero gains for unmeasured states was 0. 7494 x 10 . The 

O 

quadratic cost (J**) for this set of gains was 0. 9344 x 10 compared to 

g 

0.8664 x 10 for perfect sensing. The differences in these costs are 
sufficient to yield a meaningful test of the gradient- search method even 
though the values of J* are small. 
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Detailed equations for the gradient iterations are given in Appendix A. One 

step of the search based on minimizing the quadratic cost J** lowered this 

o _ 5 

cost to 0. 9265 x 10 . The J* cost, however, went up to 0. 2188 x 10 . 

This was because the quadratic Q(t) used for the perfect-sensing controller 

is not necessarily the best for the constrained system as indicated above. 

The cost J* can be minimized directly in a gradient search as shown in 
Appendix A. This method proved very successful in reducing the J* cost 
from 0. 7494 x 10 _6 to 0. 1266 x 10" 8 to 0. 1385 x 10" 12 to 0. 110 x 10 _l7 in 
three iterations. 

Initial conditions for the gradient iterations were obtained by approximating 
the Kalman Estimator with time-varying gains. Both methods (DC and 
averaging over frequency) were used with data at 33 instants of time, 5 
seconds apart. Linear interpolation was used between these points. The 
first method [ K(t) = Y(0, t)] proved to be quite poor because the pitch anc 
pitch rate states are of higher frequency. Thus these DC gains were too 
far off. The resulting quadratic cost was about four orders of magnitude 
higher than for perfect sensing. The probability of mission failure was 
certainty (J*=88). The second method proved to be somewhat better. The 
quadratic cost was up less than one order of magnitude. However, 
probability of mission failure was still certainty (J* = 2). It was still close' 
enough for the gradient search. One gradient step was taken, and the J* 
cost was lowered to 0. 354. 

For problems, where the unmeasured states are important, which was not the' 
case here, this method of finding initial conditions appears to be satisfactory. 
Use of more data points in time and frequency would probably give better 
results. The weighted averages were computed using only 12 values of w 
between 0 and 10 radians per second. 
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As a severe test, the gradient method was tried with the d-c gains as initial 
gains. Using appropriate normalization of the gradient vector and using the 
initial quadratic weights, we were able to reduce the J* cost to 0. 37 3. in 
three steps. For the first two steps, the gradient was normalized for each 
gain over all time separately. This had the effect of changing the shape of 
each gain versus time. If this hadn't been done, and the gradient had been 
normalized over all the gains over all time, only the shape of one gain 
(drift) would be changed, since the gradient for that gain always dominated 
the others. 

For the first two steps (J* = 88 to J* = 0. 525) the original quadratic weights 
were used to drive the response covariances down to a reasonable level. 

The probability calculations were not valid for the original covariance 
values. Thus, minimizing the J* cost was unsatisfactory because the 
gradient was invalid. The quadratic cost in the first two steps was reduced 
four orders of magnitude. 

After the first two steps, the J* cost was low enough for valid calculations 
of the partial derivatives of J* with respect to covariance responses so that 
they could be used for the quadratic weights. One more step reduced the J* 
cost to 0. 373. * 

The gradient iterations were terminated with this set of gains since previous 
iteration had been successful in reducing the cost to a desirable value from 
such a value of J*. 

Equations for simplification of the time-varying nature of the gains were 
derived and appear in Appendix B. 

The time-varying gains obtained with gradient iterations from the three 
initial conditions chosen are shown in Appendix C. These figures demonstrate 
that the second and third iterations from the perfect measurement gains were 
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of a fine-tuning nature. The figures also demonstrate improved behavior 
starting from the averaged approximation rather than the DC approximation 
to the Kalman controller. 


CONSTANT COEFFICIENT EXAMPLE 

The test vehicle considered is a 23rd-order model of the C-5A pitch dynamics. 
The model consists of the short-period mode, six flexure modes, two second- 
order Kussner lift-growth modes, two wind states, and three first-order 
actuators. The three control surfaces are inboard elevators, .ailerons and 
spoilers. The .response vector consists of stresses and stress rat.es .at five 
stations and normal accelerations at four stations. A more .detailed 
description of the C-5A model is given in Appendix E. A possibly practical- 
controller form was chosen on the basis of previous analysis pf this model 
performed on the LAMS program (Ref. 10). This form is shown in Figure 5. 



Figure 5. Controller Configuration for C-5A 
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Attempts to use the Implicit Function Method described above required 

small steps in X (AX = -0. 05). Investigation revealed that the gradient of J 

with respect to gains in K 1 (denoted by vJ) was highly sensitive to changes 

in gains. Small percentage changes in elements of K'*' gave rise to orders- 

2 

of-magnitude changes in components of vJ. The matrix y J was rather ill- 

10 

conditioned with a range of 10 for the magnitudes of its eigenvalues. The 

original algorithm w'ith AX = -0. 05 was unsuccessful at X = 0. 80. The predicted 

controller was unstable. The algorithm was modified to include damping, and 

two gains which appeared to oscillate in a diverging manner were held constant 
3 

(put in K ). With these modifications, satisfactory performance was obtained 

to X = 0. 50. The gains at X = 0. 45 were unsatisfactory. They yielded a non- 

2 

definite matrix of second partials V J and consequently unsatisfactory gains 
at X = 0. 40. 

The Implicit Function Method was discarded in favor of extrapolation followed 

1 

by gradient corrections. A predicted set of gains for K at X = 0 was com- 
puted from the equation 

K 1 (0) = 20 K 1 (.95) - 19 K 1 (l. 0) (31) 

Ten gradient corrections were made from this predicted gain. The magnitude 
of the largest component of yj for the resulting gain was less than the magnitude 
of the largest component of VJ at X = 1. 0. This was considered to.be the 
accuracy that could be expected from the numerical calculations. 

Finally, the incremental gradient method was used with a slightly different 
controller configuration. This configuration consisted of allowing flexure 
feedbacks to the elevator and spoilers. For this calculation AX was chosen 
to be -0. 10, and four conjugate gradient steps were used for correction at 
each value of X. This method proved to be quite successful for the case 
studied. 
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Numerical results are shown in Tables 1 through 4. In Table 1, the behavior 
of the Implicit Function Method algorithm for the C-5A is displayed. Without; 
damping, the algorithm failed at X = 0. 8. The modified algorithm to include 
damping was slightly better, but failed at X = 0. 70. A second attempt with 
the modified algorithm was made in which two gains corresponding to flexure- 
mode displacements were constrained to be constant. Greater success was 
achieved in that failure did not occur until X = 0. 35. However, the last step 
was very poor, and simple extrapolation proved to be much better. This 
motivated the extrapolation to X = 0 followed by successive gradient steps. 

The results appear in Table 2, where step number zero is the result of 
extrapolation from X = 0. 95 to X = 0. The gradient steps were terminated 
at the point where the magnitude of the gradient vector, as measured by its 
largest element, became less than the magnitude of the gradient at X =1. 

The gains for this controller are given for comparison with the original 
gains and those from extrapolation in Table 3. Finally in Table 4,, the per- 
formance of the Incremental Gradient Algorithm is listed. Flexure feed- 
back to all control' inputs was permitted in this case. The cost in the second 
column is that obtained before gradient correction while column three 
indicates the cost after gradient correction. 

The magnitudes of costs and gradients listed in Tables 1 through 4 may 
appear somewhat large considering the fact that the gradient for the optimal 
control should be zero. However, these magnitudes are arbitrary in the 
sense that the scaling of J and hence of VJ is arbitrary. To obtain a calibra- 
tion of J with respect to physical considerations, the value of J corresponding 
to no stress and stress rate costs and standard deviations of 0. 1 g for the 
accelerations at the four stations is 8500. Significance of the magnitude of 
vJ computed for the optimal controller with complete measurement feedback 
may be deduced as follows. Suppose that K* denotes the optimal gain matrix 
which will be considered as a gain vector. In a sufficiently small neighborhood 
of K* the cost may be represented as 



Then vj J K = v 2 J | K; ,.(K-K-). Now suppose K-K* = EK* 

where E is a diagonal matrix. The diagonal elements of E may be considered 
as relative errors in the components of K* and should be small if K approxi- 
mates K*. For the computed optimal controller each element of E was less 
than 0. 019 in magnitude indicating an accuracy of approximately 2% or better 
for each gain. 

Table 1. Behavior of Implicit Function Method Algorithm for the C-5A 


Original Algorithm 

Modified Algorithm 




Unconstrained K's 

Constrained K' s 

X 

J 

max | (vj ). 1 

n 

maxj(vJ) i j 

J 

max|(vJ).. | 

1.00 

14377 

.28 • 10 5 

14377 

. 28 • 10 5 

14377 

. 28 • 10 5 

.95 

14429 

. 88 • 10 5 

14429 

. 43 * 10 5 

14427 

. 14 • 10 5 

.90 

14595 

. 98 * 10 5 

14590 

. 45 • 10 5 

14582 

.41 • 10 5 

.85 

14884 

. 42 • 10 6 

14862 

.91 • 10 5 

14848 

.28 • 10 5 

.80 



15614 

.23 • 10 6 

15233 

. 35 • 10 5 

.75 



16092 

. 34 * 10 6 

15749 

. 12 • 10 5 

.70 





16416 

. 31 • 10 5 

.65 





17248 

. 57 • 10 5 

.60 





18275 

. 53 • 10 5 

. 55 





19470 

. 52 • 10 5 

. 50 





21170 

. 18 • 10 6 

.45 





24366 

. 43 • 10 6 

.40 





86984 

. 95 • 10 7 

. 40 a 





26482 

. 53 • 10 6 


indicates extrapolated set of gains 
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Table 2. Extrapolation plus Gradient Algorithm 


Step 

Number 

J 

10 5 * max j( vJ)-| 

0 

97125 

. 51 

1 

94789 

34 

2 

93389 

18 

3 

92373 

22 

4 

88926 

13 

5 

88644 

10 

6 

88445 

6. 3 

7 

83880 

0. 5 

8 

72191 

3. 8 

9 

70685 

0. 9 

10 

70679 

0. 275 


Table 3. Gain Comparison 


Gains 

Controller 

Optimal with 
Complete Feedback 

Extrapolated 
(step 0) 

Final 
(step 10) 

K i.i 

- ,.02983 

- .00780 

- .00545 

- K l,2 

- .01574 

. 0103 

. 0114 

K l,3 ' 

- .2388 

.00520 

. 00899 

K l,4 

.06796 

.0900 

. 0842 

K l,5 

2.837 

-1. 06 

-1. 06 

K l,6 

. 3137 

.048 

. 100 

K l,7 

-1. 396 

-4. 90 

-4. 90 

K 1,S 

- .4093 

- . 175 

- .023 
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Table 3. Gain Comparison (Continued) 


Gains 

Controller 

Optimal with 
Complete Feedback 

Extrapolated 
(step 0) 

Final 
(step 10) 

K l, 10 

. 685 

. 345 

. 005 

K l* 12 

.07228 

.0723 

. 0720 

K l, 13 

. 9918 

.491 

.491 

K l, 14 

- .6193 

- .485 

- . 160 

K 2,l 

. 008698 

. 00250 

.00158 

K 2* 2 

. 004714 

- .00850 

- .00110 

K 3,l 

. 001872 

- .00073 

. 00013 

K 3* 2 

.001119 

.00072 

. 00633 


Table 4. Performance of Incremental Gradient Algorithm 9, 


X 

J 

(predicted gains') : - 

J 

(corrected gains) 

1 . 0 

14377 

14377 

. 9 

15380 

15186 

.8 

17375 

16387 

.7 

18664 

■17479 

. 6 

19394 

18870 

. 5 

25517 

21308 

.4 

24179 

23181 

.3 

26086 

24832 

.2 

27202 

25867 

. 1 

28146 

27194 

0 

30611 

28015 


a J for the free aircraft is 116* 193. 
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SECTION in 

PARAMETER VARIATION CONSIDERATIONS 


The goal of implicitly including parameter- variation constraints within 
quadratic optimization formulations motivated study in two directions. One 
phase of the study was aimed at discovering properties of the optimal per- 
formance surface over a segment of parameter space and utilizing such 
properties to find optimal insensitive controllers. The other phase of study 
was aimed at reducing sensitivity to parameter variations by introducing 
proper compensators into the controller. Promising results of both 
theoretical and practical value were obtained in the first phase. Although 
intuitively appealing, the problem formulation of the second phase appears 
to be of questionable practical value. 


PROPERTIES OF THE OPTIMAL PERFORMANCE SURFACE 

Consider the scalar equation x = fx + gu + q with performance functional 
J = E J (4x +u )dt where q is a unity white-noise input. The coefficients 
of the System, f and g, are constants with values in the rectangle R defined 
by 0 ^ f ^ 2 and 1 ^ g ^ 5. For a controller, u = kx, the performance 
functional is a function of the three parameters, k, f and g; i. e. , J = J(k, f, g). 
The optimal control for any (f, g) in R is u = k*(f, g)x where k*(f, g) minimizes 

J(k, f, g). For this simple example, the solution for k*(f, g) is easily determined 

V o T 

f + 4g ]g . In this case J(k*(f, g), f, g) can be obtained 
explicitly as -k*(f, g)g - ^. This surface over the rectangle R, as depicted in ' 
Figure 6, has a maximum over the corner (2, 1) and a minimum over the 
opposite corner (0, 5). Further calculations show. that the surface J(k*(l, 3)f, g) 
corresponding to performance with the optimal controller for the midpoint 
of R highly accentuates the peak at (2, 1) while providing nearly optimal' 

PEECEDEto pa® blamr 
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performance for points close to the midpoint of R as shown in Figure 7. This 
controller appears; to be sensitive to parameter variations. To alleviate this 
sensitivity,, a controller could’ be chosen to minimize E f (J(k, f, g)} assuming 
a uniform probability distribution for the parameters f and g in the rectangle 
R. With this simple example, this can be accomplished yielding a value of 
k = -2. 82.. Indeed the peak of J(k, f, g)- which occurs at (2,.l) is lower, as 
shown in Figure 8. But the- controller which minimizes the maximum of 
J(kj f, g ) over R is the optimal controller for the point (2, 1). Any other k 
will give a higher value of J(k, 2, 1). Furthermore J(k*(2, l)f, g) is less than 
or equal to- J(k*(2, 1), 2, 1) as shown in Figure 9. This controller is the 
least, sensitive to parameter variations defined by R. 

The following two, theorems indicate that the phenomenon noted in the above 
example is not mere happenstance. The first theorem is an existence 
theorem providing a sufficiency condition for optimal insensitive controllers 
locally, that is, for sufficiently small variations. This theorem also implies 
that such controllers are optimal controllers corresponding to boundary • 
points of admissible parameter variation sets. The second th eorem states 
a necessary condition for an optimal controller corresponding to a point in 
the boundary of the domain 1 of admissible parameter variations to be an 
■optimal insensitive controller. Here, the expression optimal insensitive 
controller refers to a controller which is optimal for some admissible value , 
of the parameters which minimizes the maximum of the -performance over 
'the range- of admissible parameter values. 

• Theorem 1— Consider the system x = F(p)x + .G(p)u, x(0) = x_ and 

f»00 T ^ 

an associated performance functional J(u, p) = J! (Hx+Du) Q(Hx+Du)dt 
where- p is a- vector of : parameters.. Let J*(p) = mm J(u, p). Suppose 
p Q is ,a point with the property that VpJ*(p Q ) t 0* Then there 
exists an-e >0- such that the control u*(p ) which minimizes 
J(u, p Q ) -also^minimizes the maximum of J(u, p) with respect to 
p in an e-ball with p Q on the shell (boundary) of the -ball. 
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Figure 8. Performance Surface for Optimized 
Expected Parameters 


Figure 9. Performance Surface for Optimal 
Insensitive Controller 






Proof: For any e > 0 let B_(p ) denote the e-ball with center at 

GO 

p 0 - e y* ( Po>- i - e -> 

B e = jp : P = P 0 -eV p J*(p 0 )+eri, in | £ |V p J*(p 0 >| 


Also define M(u; p , e) to be max J(u, p). Then 


o 


P eB e (P 0 ) 


M(u; p^,e) s J(u, p o ) s J*(p n ). For peB g (p Q ) we may express J(u*(p n ), p) 


o 


as 


J(u*(p Q ), p) = J(u*(p o ), p o ) + ^ J(u*(p o ), p) 


. (p-p Q ) + 

P=P„ 


l/2(p-p o ) T 


5 2 J<u*(p o ), p) 



(P~P 0 )+o( e ) 2 

P=P G 


= J*(p Q ) + g* (p-P Q ) + 1 /2(p-p o ) T H(p-p o )-+ o(e 2 ) } H> 0 

(32) 


Note that g = V J(u*(p ),p) 
P ® 


= V J*(p) 

p 


P=P, 


P = P, 


For ps B. and e sufficiently small the only possibilities for 
extreme points of J(u#(p ), p) are 

(1) approximately p ~H if this point lies within or 

(2) points on the shell of B . 
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The point near p o -H -1 g is a minimizing point. Therefore 
maximizing points lie on the boundary of B e . 

The problem of extremizing J(u*(p 0 ), p) subject to p = P 0 -e g +eT l 
with | r| | = | g | may be treated with a Lagrange multiplier as 
minimizing 

H = J q + g * (p-p Q ) + l/2(p-p o ) T H(p-p o ) + o(e 2 ) + \(rr r\~ g. g). 

2 2 i 

, This yields 0 = e g +Ah-e H(g-rj) + 0(e ) and | r| | = | g | as 
necessary conditions. For e small these conditions imply that 

X = ±e [l + o(l)] and n = + g[l + o(l)]. 

The bottom signs yield the maximum and the top signs yield the 
minimum. The exact solution for the bottom signs is X = -e, T|= g 
which describes the point p . Thus on B„, J(u*(p ), p) <y(u*(p ),p ) 
= J*(p ). Hence 

M(u*(p Q ); p Q , e ) = J*(p 0 ) ^ M(u; PqJ e) 
which was- to be proved. 

Theorem 2 - Consider J(u*(p Q ); p) = J*(p. ) + g • (-p-p) + o( |p-p |) 

with g = V J*(p) 1 . Let Q denote a closed convex set with non- 

s P 'P = Po 

empty interior in the parameter space. Suppose p is a point in 
the boundary of Qwith the property that 

J(u*(p Q ), p Q ) = max J(u*<p o ), p). 
pefi 

Then g must be an external normal to Q at p Q . 
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Proof: Assume that g is not an external normal to Oat p 

o’ 

l. e. , there exists a p 1 in Qsuch that (p„ -p ) • g = c >0 Since 

J- 1 o ° 1 

0 is convex p(A) = P Q + ^(p 1 -P 0 ) lies in Q for 0 ^X £1 and 
(p(\)-p o ) • g = X (P^~P Q ) * g = Xc^. Thus 

J( u *(p Q ), p(\)) = J*(p o ) + (p(X)-p o )* g+ o( I p{X)-p o |) 

= J"<P o ) + XCj + o(X) > J*(p q ) (33) 

for X sufficiently small. This contradicts the hypothesis that p 

has the property that J(u*(p ), p ) = max J(u*(p ), p). 

° pe n ° 


To utilize information about performance surfaces in determining insensitive 
controllers, approximations and computational techniques are required. For 
this purpose, consider the system 

x = F(p)x + G(p)u, x(0) = x q 

and a performance index 
00 

J = J <Hx+Du) T Q(Hx+Du)dt 
o 

Suppose J*(p) denotes the optimal performance surface, J(K*(p), p), 
corresponding to u = K*(p)x. For an arbitrary controller u = Kx, let J(K, p) 
denote the corresponding cost surface. Define the "error" index as the 
difference 

e(K, p) A J{K, p) - J*(p) 


39 



This error function may be expressed as an implicit function of K and p as 
e(K, p) = Tr jAX Q [ - Tr jpX Q J (34) 

T 

where X = x x £ A satisfies 
o o o * 

(F+GK) T A + A(F+GH) + (H+DK) T Q(H+DK) = 0 
and P satisfies 

F T P + PF + H T QH - PG(D T QD)" i G T P - PG(D T QD)'' 1 D T QH 
- H T QD(D T QD) _:L G T P - H T QD(D T QD)" 1 D T QH = 0 

A is an implicit function of K and p and P is an implicit function of only p. 

So J(K) and J*(p) can be approximated by truncated Taylor series expansions. 
It is assumed that the parameter variations of F and G are moderate, say 
10 or 20 percent variation from the nominal value. Therefore, we might 
not need higher-order partials of A and P with respect to K and p. At most, 
it is assumed sufficient to expand J and J* in Taylor series up to the second- 
order partials of A and P. It might turn out for a practical system that the 
first-order expansion of J and J* would be sufficiently accurate. In any case 
the error function e(K, p) can be explicitly expressed in terms of K and p in 
a reasonably accurate form.. Its derivation is given in Appendix D. A 
simple- geometrical interpretation of the approximations of J and J* is shown 
in Figure 10, where p is a one- dimensional parameter, and K is fixed. 

In Appendix D the following facts are established. 

(1) e(K*(p), p) = 0 

This implies that the surfaces J*(p) and J(K*(p), p) touch at p. 
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J*(p) = TRUNCATED TAYLOR SERIES APPROXIMATION OF J*(p> 

Figure 10. A Geometrical Interpretation of Taylor Series 
Approximations of J(K,p) and J*(p) 

(2) V e(K*(p), p) = 0 

P 

This implies the surfaces are tangent at p. 

(3) V K e(K*(p), P) = 0 

This is the statement that K*(p) satisfies the first-order condition 
of optimality. 

( 4 > v KK e(K*(p), p) > 0 

This, together with (3) states that K*(p) is optimal. 


These properties along with computational techniques described in Appendix D 
for obtaining truncated T.aylor Series approximations make it reasonable to 
compute optimal insensitive controllers for realistic systems. 


Results similar to those above were derived by Salmon (Ref. 20). 
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As a test vehicle for this approach, the C-5A model was modified to include 
only rigid-body and three flexure degrees of freedom and first-order 
actuator dynamics. Variations in the natural frequencies of the flexure 
modes were introduced as the parameter variations. Numerical results 
are presented in Appendix E. 


INSENSITIVITY VIA COMPENSATORS 

This phase of the study was aimed at relating the structure of the con- 
troller to the sensitivity of the controlled system. For a system represented 
by the familiar equations 

x = Fx + G^u + Ggtl 
r = Hx + Du 

and the associated performance functional 
J = E {Jr" 1, Qr dt } 

where the matrices F, G^ G 2 , H, D and Q may depend on a vector of 
parameters, p, it is known that the optimal controller is given by u - Kx 
where K may depend on p . The equation u = Kx has a simple geometric 
interpretation. The equation defines a .subspace of the (x, u) space. With 
x and u scalars this interpretation is reduced to: u = Kx defines a line in 
the (x, u) plane. Then optimality would be equivalent to motion in the (x, u) 
plane being constrained to the proper line, u = Kx, and (for the regulator 
problem) that undisturbed motion be directed toward the origin. Thus for 
a simple scalar system optimality might be- depicted as in Figure 11 where 
the "optimal" lines for three values of p are sketched. 
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Figure 11. Optimal Motion in the' (x, u) Plane 

Compensation of the feedback may be represented in differential equation 
form by augmenting the equations with 

u = Kx + K q V 

v = w 

where w is treated as the control input and the state vector is augmented to 
include v. For the augmented system, optimality may be approximated by 
asking that the motion of the augmented system by approximately in the 
"optimal" subspace. Then a performance criterion could take the following 
form: motion originating in the "optimal" subspace should 'remain in that 
subspace and be appropriately stable, motion originating outside the "optimal" 
subspace should approach the "optimal" subspace in an appropriate manner. 

In differential-equation terminology, this criterion can be expressed simply 
as requiring the "optimal" subspace to be a stable invariant manifold. 
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The^ingle scalar system x = Fx + Gu, x(0) given, with performance index 
J = J (Qx +u )dt was analyzed to gain insight into the effectiveness of com- 
pensator structure in alleviating sensitivity to parameter variations. For 
numerical calculations it was assumed that F = l+2p, G = 3+4p; Q = 2. 2 5 
and that -0. 6 £ p s + Q. 6. 

The optimal control for any fixed value of p is u v K#(p)x with 

k* = -Cf +fc 2 +qg 2 3g" 1 ' 

The desired behavior of the controlled system was assumed to be that a 
closed-loop system principal axis match the line u = K*(p)x to a high degree. 
The Laplace transform of the system with w = K^x + K 2 v is (without loss of 
generality K q is set equal to 1),: 

(s-F)X(s) = x(0) + G[KX(s)+V(s)] 

(s-K 2 )V(s) = v{0) +K 1 X(s) 

U(s) = KX(s) + V (s) (35) 

From these equations 

v(s) - x(0)Kj + v(0)(s-F-GK) 

= K + 

X(s) x(0)(s-K 2 ) +Gv(0) 

The requirement that u(t) = K*x(t) if u(0) = K*x(Q) yields-the equation: 

+ (K 2 -F-GK)(K*-K) -*G(K*-K) 2 = 0 

The left-hand side may be expanded as a Ta-ylor' Series In p -about p = 0 
yielding 
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K x + (K 2 -1-3K)[K*{0 )-k3 - 3[k*(0)-K] 2 

+ p {- (2+4K>[K*(0)-K] - 4 [k*(0)-kJ 2 + (K 2 -l-3K-6[K*(0)-K])[K*(0)] , f 

+1/2p 2 {(K 2 -1-3K-6[k*(0)-K])[K*(0)]" - 6(Ck*(0)]') 2 

- (2+4K+8Ck*{0)-k 3)[K*(0)]'} + higher order terms (36) 

The control parameters K, Ki and K 2 can be chosen so that the constant term 
and the coefficients of p and p2 are zero to yield desired matching of the 
principal axis to u = K*(p)x for p near zero. Stability requirements of the 
resulting closed -loop system are not guaranteed. This may impose further 
constraints on the control parameters and, in so doing, reduce the degree 
of approximation. 

For the numerical example it was found that matching to third-order terms 
in p could be achieved with K = -2. 08, = -1. 03 and Kg = 0. 29. The 

corresponding closed-loop poles for p = 0 are -4. 61 and -0. 34. The pole 
at -4. 61 is the desired value for the optimally controlled scalar system 
with p = 0. Thus the compensator pole now is the dominant pole, and the 
corresponding principal coordinate displays the dominant dynamics. This 
is undesirable. The principal axis which approximates u = K*(p)x is an 
unstable axis. Motion originating near this axis does not remain close to 
this axis. Figure 12 is a sketch showing this kind of behavior. 

One other case was considered in which K was set equal to zero. The 
remaining control parameters and Kg were chosen to eliminate the 
constant and first-order terms in p. The value of the gains were = -8. 15 
and Kg = -48.2. The corresponding closed loop poles for p = 0 were -4. 61 
and -42. 7. In this case the axis u = K*x was stable and the resulting con- 
troller exhibited desired properties.,. Characteristics of this controller 
are given in Table 5. The final column of this table characterizes the 
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Figure 12. Phase-Plane Portrait Indicating Unstable Subspace u = K*x 

behavior of the optimal controller based on p = 0 for comparison. Columns 
2 and 3 indicate good correlation between the dominant principal axis and the 
line u = K*<p)x. Columns 4, 5 and 8 show eigenvalues with good matching 
between columns 4 and 5. On this basis the compensated system appears 
to be less sensitive than an uncompensated system. However, the improve- 
ment may not be worth the added complexity. 

Thus, although this method of alleviating sensitivity has some appeal, its 
utility would have to be determined by application to a more realistic and 
more complicated example. Rather than pursue this, it was considered 
more advantageous to study the performance surface for a realistic example. 
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Table 5. Optimal Control and Principal Coordinate Dependence on p 


1 


3 

4 

5 

6 

7 

8 

p 

mm 

V 

F + GK* 

P 

-v 

X 

F + GK*(0) 



-2. 08 

-8. 39 

-9. 03 


-37. 0 

-7. 90 



-2. 00 

1 

-a 

• 

i »iA 

CO 

-7. 38 


-39. 1 

o 

CO 

• 

CD 

1 


-1. 91 

-1. 93 

-5. 87 

-5. 92 


O 

T— l 
1 

-5. 71 


-1. 89 

-1. 90 

-5. 23 

-5.25 

-12. 7 

l 

h-* 

• 

CO 

-5. 16 

0 

-1. 87 

-1. 87 

-4.61 

-4.61 

-14. 6 

-42. 7 

-4.61 

i 

-1. 84 

-1. 84 

-3. 98 

-4. 00 

-17. 0 

i 

CO 

• 

CD 

-4.06 


-1. 80 

-1. 82 

-3. 36 

-3. 40 

-20.4 

-44.3 

-3.51 


-1. 65 

-1. 78 

-2. 11 

-2.29 

-32. 8 

-45. 8 

-2.42 


-1. 20 

-1. 74 

- . 92 

-1.24 

-78. 1 

-47.2 

-1.32 


^Principal coordinates are = Xg + and z ^ = x^ + q^x 
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SECTION IV 
SENSOR CHOICE 


An important problem in the design of practical controllers for large flexible 
vehicles is the choice of types and locations of sensors to generate the feed- 
back signals. The problem arises because it is economically prohibitive in 
most applications to sense all system states, particularly the rates and displace- 
ments of numerous flexure modes. Moreover, experience on previous 
flexure control designs has shown that adequate performance can be obtained 
with significantly fewer sensors than state variables, suggesting that 
economics vs. performance comparisons favor smaller sensor complements. 

In principle, the sensor-choice problem is solved by selecting a complement 
of instruments which exhibits the most desirable cost/performance tradeoff, 
assuming, of course, that the instruments are utilized in an optimum 
fashion, i. e., by locating them optimally along the vehicle and by using 
optimally the information thus provided. It is in these areas that present 
methodology fails. We have too little computing power, or too cumbersome 
computing methods, to determine optimal- locations of sensors and to utilize 
their information optimally with a controller subject to various simplicity 
constraints. Cost/performance comparison for several competing sensor 
complements are thus difficult to come by. 

The research reported in this section is intended to improve this situation. 

Its general objectives are twofold: (1) to increase understanding of the 
relationships between sensor complements and the performance capabilities 
they offer and (2) to improve computational methods for locating and 
utilizing the instruments. 


PRECEDING PAGE BLANK NOT FILMED 
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MATHEMATICS OF THE SENSOR-CHOICE PROBLEM 


As usual, we are given a mathematical model of the following form: 

• System: 

x = F(t)x + GjttJu + G 2 (t)n (37) 

• Responses: 

r = H(t)x + D(t)u (38) 

• Performance Measure: 

T rp 1 

J ■= E j J* r(t) Qr(t) dtf (39) 

*o 

where all symbols -correspond to -previous definitions. Then we are given 
a collection of measurement instruments (rate and acceleration sensors, 
stress sensors, air data sensors, etc. ) from which we may select any 
number and combination to mount anywhere on the airframe. The outputs 
of these instruments are to be used for control. That is, let Q denote a 
particular subset of instruments, and let vector y denote the locations of 

J.T_ 

each instrument (i. e., y^ is the location of the i sensor). Each location 
can be represented by a scalar number if fuselage and wing positions are laid 
end to end. Then the instrument outputs, at least for linearized models, will 
be 


z ■= M(t, y, Q)x 

and the controller must be of the form 


(40) 


u(t) = L(t, z) 


(41) 



In these latter equations, M(t, y, q) denotes a matrix of dimension (mxn), where 
m is the number of sensors in the set Q, and the function L(t, z) is a linear con- 
trol law which is usually constrained to exhibit certain attributes of simplicity. 
For instance, the nondynamic form, L(t, z) = K(t)z, or a dynamic form with 
low-order compensation are desirable. 

With measurements (40) and controller (41) substituted into equation (37), the 
performance measure (39) becomes a function of n> y and L, i. e., 

J = J(n,y, L) (42) 

This expression shows explicitly why cost/performance comparisons of sensor 
complements are difficult to get. To evaluate the performance capability, 

<T'(n), of the set of instruments q (whose cost at a given level of reliability is 
presumably known), we must specify both the sensor locations and the control 
law which is best for Q, i. e. , 

J^Mfi) = min min J(Q, y, L) 

y l 

= J(a y*(n), l-(.,.,q)) ' (43) 

With current capabilities, the simultaneous optimizations of y and L, with L 
subject to simplicity constraints, are expensive indeed. Just to get a feel for 
the magnitude of the problem, suppose we treat a single fixed-time point of 
the general nonstationary situation. Suppose further that a nondynamic con- 
troller of the form L(z) = Kz is sought. Then we could use one of the parametric 
methods in Section II to compute optimal controllers, L, as functions of y and n 
and we could imbed these computations inside an iterative Newton-Raphson or 
gradient loop to optimize y. Diagrammatically, the procedure would like like 
Figure 13. 

Depending upon the order and complexity of the problem, the inner loop of this 
algorithm [block (1)] consumes anywhere from one-half to three hours of 
computing, say on an SDS 9300 machine, t The outer loop updating step 

-(■The one -half-hour figure was obtained on a 20th-order single flight con- 
dition optimization for the F4 aircraft (Ref. 6), while the three-hour 
figure was obtained for the 23rd-order C-5A example discussed in Section 
II. 
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Figure 13. Basic Sensor-Choice Algorithm 

(block (2)] can be performed either by repeated executions of block (1) at 
perturbed values of y or by the solution of coupled state covariance and 
costate equations which result from analytic evaluations of the first and 
second variations of J. For first partials alone, the latter equations com- 
prise a system of n{n-l) linear equations. The computations in both blocks 
are thus extremely time consuming — and, of course, they must be repeated 
for as many outer-loop iterations as are required for convergence. 

While it is thus apparent that the algorithm of Figure 13 is ill-suited to the 
general sensor-choice problem, it should be noted that it has been seriously 
proposed for more specialized versions of the problem (Ref. 11) and has 
been successfully run on a specialized version of the "dual" problem — force 
producer choice (Ref. 5). In each case, the specializations were such as 
to greatly reduce the computational loads of blocks (1) and (2).. Namely, it 
was assumed that controller JL, is subject to no practicality constraints. 

This leads in block (1) to the controller L*(q, y) which consists of a cascaded 
Kalman filter and Kalman controller (a system with n 11 order compensation). 







Both the filter and controller are comparatively easy to compute. Moreover, 
the variations in block (2) also turn out to be simpler. This is because the 
Kalman controller is independent of y in the case of sensor choice, and the 
Kalman filter is independent of y in the case of force producer choice. The 
previously coupled covariance and costate equations uncouple and leave 
comparatively simple equations for the variations of J. 

i 

Since we are concerned here specifically with constrained controllers, the 
simplifications offered by these special cases are inapplicable and the 
algorithm of Figure 13 offers little hope toward solving the sensor-choice 
problem efficiently. The research reported here has therefore been directed 
along different lines. Instead of relying on the quadratic cost, J*(Q), as a 
performance measure of a sensor complement, the research was aimed toward 
development of alternate quality measures with two key properties: 

(1) They should be easy to evaluate, for use in computational 
algorithms such as Figure 13. 

(2) They should provide meaningful indications of performance 
capabilities offered by a set of sensors. Ideally, they should 
exhibit strong correlation with the cost J*(q). 

This line of attack follows the approach taken in Reference 3, where two 
alternate quality measures were proposed, though neither proved wholly 
successful in satisfying property (2). 

Like Reference 3, the present research has no final wholly successful 
answers to report. Two promising quality measures were developed, both 
based on the pole -placement capabilities of the sensor complement. As 
such, they strictly apply .only to stationary problems, but can be used on 
nonstationary ones via frozen-point procedures. Most of the development 
work has been concerned with establishing theoretical equations and 
verifying the computational characteristics required by property (1). 
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Results of these efforts are described below. They are general in nature 
and find .application outside of the sensor choice problem in areas such as 
multiloop controller design, initiation of computational algorithms, and 
sensitivity analysis. As far as property (2) is concerned, several basic 
performance capabilities (e.g., stability, modal characteristics) follow 
directly from the measures. Correlations with J*(Q), however, and the 
practical utility of the measures remain to be established on realistic 
design problems. ■ 


POLE -PLACEMENT QUALITY MEASURES 

The basic notion pursued here is to examine the freedom offered by a sensor 
complement in locating closed-loop eigenvalues as a possible source of 
quality measures. This is motivated by three considerations. First, 

■ classical experience with root loci and frequency domain design techniques 
.provides tested insightful relationships between the performance capabilities 
of a controlled system and the closed-loop pole arrangements permitted by 
sensors. Such notions as stability, frequencies- of oscillation, damping of' 
.individual modes of response and dominance are all apparent from the pole 
constellation. Second, there is a fundamental connection between pole 
placement and the concept of controllability. With full state feedback, 
•freedom- to assign all n system poles arbitrarily has been shown mathe- 
matically equivalent to the condition of complete controllability (Ref. 12, 13). 
'Hence it appears fruitful to investigate any restrictions of this pole 
assignability property which are imposed by partial state feedback. Third," 
in certain special control problems, . namely single input, the state dependent 
terms of the quadratic cost (39) can be expressed uniquely in terms of 
closed loop poles, with no dependence on gains or measurements used to* 
bring them about. Though this property' fails for multi-input systems and 
ignores control- dependent terms, it provides a potential- link between 
quality measures based on pole placement and the measure J*(q). 
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As shown in References 14 and 15, the primary effect of limited 
state measurement is to reduce the total number of poles which can be, 
assigned arbitrarily — from the maximum of n poles' with full state measure- 
ment to a maximum number, p , which will be defined later. The effect 

in 3.x 

does not confine each pole to some subregion of the complex plane, as one 
might expect from classical root locus plots. Rather, P max poles may be 
placed anywhere in the plane (subject to certain nonsingularity constraints), 
with the remaining ( n “P max ) poles implicitly determined. These observations 
hold for direct feedback of the measurement with a nondynamic controller 
= Kz. If compensation is allowed, the number of arbitrary poles can be 
increased by an amount also defined later. 

These properties suggest the following quality measures for sensor comple- 
ments: 

• The number of poles, P max > which can be placed arbitrarily 

• Some measure of deviation (say quadratic) of the unassignable 
poles from specified desirable positions, given that the 
assignable poles are placed in desired spots 

While these are only verbal definitions of what must eventually be analytic 
measures, their potential is apparent. The first serves as a gross differenti- 
ator between competing sensor complements. If performance requires that 
p poles be positioned accurately, then all combinations of instruments which 
fail to have this capability can be eliminated immediately. Computations 
should be simple, since the measure depends only weakly on sensor locations, 
y. The second measure provides finer differentiation between the remaining 
complements. It evaluates the amount of pole placement deterioration in 
noncritical modes of response, given that the critical modes are adequately 
controlled. This measure should depend strongly on y and could thus be 
used to locate the sensors by means of an algorithm such as Figure 13. For 
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both measures, the pole constellations generated by quadratic -optimal un- 
constrained controllers could be used to define critical and noncritical modes 
of response and desirable pole locations. 

All this is conditioned, of course, on computational feasibility — whether it 
is possible to determine the maximum number of arbitrary poles and to 
design the controller which achieves the desired critical closed-loop locations 
with computational efforts significantly below the J*(fi) requirements. Pole- 
placement equations and algorithms in the literature suggest an affirmative 
answer (Ref. 14, 15, 16, 17). However, no published procedures were 
found which completely solve the specific problems posed here. For a non- 
dynamic multi-input multi-output controller, for example, available results 

establish that p = m, where m is the rank of the measurement matrix, 
^max 

M(q, y). This number is based on analytical methods which’ reduce the 
system to single-input form before assessing pole assignability properties. 

In this study it has been shown that significantly more poles can actually be 
placed when the additional degrees of freedom offered by multiple inputs are- 
utilized. Similarly, no figure seems available for the number of arbitrary 
poles added by a compensator of specified order, and in particular, no 
procedure was found to specify the minimum compensator order required to’ 
achieve arbitrary placement of all system poles and to compute the param- 
eters of this compensator. These questions were answered in the course of 
this research. The answers take the form of a general set of" pole -placement 
equations and a computational algorithm to solve them. Both are discussed 
in the remainder of this section. The algorithm provides a computationally 
reasonable way to evaluate the quality measures proposed above. This is 
demonstrated in Appendix G with a few trial computations for'F-4 and C-5A 
flight control design problems. 
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GENERAL POLE -PLACEMENT EQUATIONS 


The poles of a closed-loop stationary system (37) with nondynamic controller 
L = Kz are roots of the following characteristic polynomial: 


min (r, m) 
P{s) = D(s) + s 
k=l 



} l m ** J k 


r. 


N 


i k a v • • £ k i 




(s) 


T 

P 
L i 


(±) K. # K 


3i^i 12*2 



(44) 


Here P(s) is the closed-loop polynomial, D(s) is the open loop polynomial of 
(37) with 11 = 0 , and the expressions N 3 ^ \ are various polynominals 

associated with the feedback loops. The summations indexed by j 2 , ... j k 


and 


l are carried out over all naturally ordered groups of k 

K 


out of r controls and k out of m measurements, respectively, and the sum- 
mation indexed by p^is carried out over all permutations of the sequence 
l l . . . , with algebraic sign taken positive for even permutations 
and negative for .odd. 


Equation (44) is a polynomial in the complex variable s and in the mr gain 
variables K, with order n in s and min(m, r) in K. Its detailed derivation is 
left to Appendix F. Here we will present only interpretations of terms and 
some procedures for computing them. We will then use shorthand versions 
of the equation to derive necessary and sufficient conditions for placement of 
p poles and to develop explicit expressions for P max * These will be shown to 
apply to dynamic as well as nondynamic controllers; they will be illustrated 
with a small example pole -placement problem; and finally, they will be used 
to devise a computerized pole -placement algorithm. 
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Interpretations and Definitions 
3l • 3k 

The terms, N „ (s) may be treated as generalized numerator .poly- 

’ ** - ** * 

nominals of the stationary multi -input system (37) with .multiple stationary 
outputs (40). For example, the collection of first -order terms (k=l) are 
the familiar numerators of the transfer -function matrix 

H(s) = -M(a.y) (si-F)" 1 G 1 ; (45.) 

Each term can therefore be computed b.y replacing , certain rows .of the matrix 
(sI-F) with rows of the measurement matrix M(q, y) and evaluating the resulting 
determinant. For the termN^, this means replacement of r the i(j) row of 
(sI-F), can it (sI-F)* 1 ^, by -(Gj).^^ times the ./ h row, M U) , ( of matrix,]^. 
This procedure assumes that each column of.matri'x G 1 has a single nonzero 
entry in positions i(j), j = 1, 2, . . . , r; which will. always b_e ; the case.for 
systems with actuator dynamics- and repr.esents.np loss, of generality^ even for 
.other cases. 


The collection .of second-order terms (k=-2) are .so-called '.'coupling numerators,." 
of the system- (Ref. 1.8) which ar.e present, whenever two.- feedback paths exist 


simultaneously. Each of these terms is computed-.by. replacing, two,.rows pf 

■ ~ l- ‘ j j-j2 ’ * * r< ' 

(sI-F) with rows of M and computing determinants. For , the replace- 
ment schedule would be (sI-f/^ 1^ replaced by r (Gj.,. and 

(sI-F) (i(; te>) replaced by -(G^^M^). * 1 


Analogous replacement schedules- apply for higher.-order numerator terms, 
with k row replacements for each k^-order term. No.. simple interpretation 
of. these terms, however, seems-to be. possible except, to say that .they arise 
whenever k feedback paths exist, simultaneously. 
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Conditions for Arbitrary Pole Placement 


For ease of manipulation, it is convenient to express Equation (44) in the 
following shorthand form: 

P(s) = D(s) + E N.(s) 0 . <K) (46) 

i 

where the various generalized numerators are designated simply by bL(s) 
and the gain functions which multiply them are denoted by 0 .(K). It is now 
our objective to determine the maximum number of roots of the polynominal 
P(s) which can be assigned arbitrarily by proper choice of K and to determine 
the fate of the remaining roots. To do this, assume that the number of 
arbitrary poles is p, whose maximum is as yet unknown. Assign these poles 
to be roots of a specified polynomial \(s), and let the remaining poles be roots' 
of the polynomial 6 (s). Then (46) becomes 

P(s) = \(s) fi(s) = D(s) + 2 N.(s) 0 . (K) (47) 

i 

and in coefficient-vector-formt we get 

A 6 + X = D + s N. 0 . (K) 
i 

0 = D - \ - AS + S N 0. (K) (48) 

i 1 

where A is a matrix and \ a vector defined as follows: 


fin coefficient -vector -form, a polynomial P(s) - p. + P 2 S + P 3 s + 
p^s +s n is represented by the n-vector P = (p^^ P 2 • • *P n ) • 
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( 49 ) 


Equations (48) represent a set of n nonlinear equations of the form f(x, 6, K) =' 6' 
which define K and 6 as functions of the specified polynomial i. e.y K = K(X)’ 
and 8 = 6(x)« According to the implicit function theorem (Kef. 8), these 
functions exist in a neighborhood of a point (X Q , S Q > K q ) if; and only if the 
rank of the "Jacobian matrix" 


♦ <V 8 o- K o> = 

df 

T 

36 

1 

1— T 
1 2K T 
i J 

x 6 K 
A o o o 


— 

1 if N. 
! ‘ 1 

d0 i 

— 


' A o 

0K t 

K o- 


<5d)t 


tThis equation makes sense in vector -matrix notation only if K is written out 
as a vector. This is assumed. 
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equals n, and they are unique if in addition the matrix is square. Existence 
of solutions K(\), §(x) are, of course, the very conditions which must be 
proved to establish arbitrary placement of the p poles, \(s). We have thus 
arrived at the following basic result: 

i 

* Theorem 1 — Locally arbitrary placement of p poles is possible 
in a neighborhood of (\ , 6 , K ) if, and only if, 
rank l>(\ o , 6 q , K }] = n 


This theorem can be used to determine p and also to examine pole place- 

max r 

ment conditions in a special neighborhood of interest (*■ , 6 , 0). To determine 

o o 

P max ' note that the rank of tjr is bounded by the inequalities 


rank 


20i - 


?0i ~ 

S N - rp 

£ rank [ijr ] s rank 

S N, rp 

i 1 dK l 

i ‘3K T 


+ rank 


(51) 


with the right-hand equality realized when all columns of A are independent 


of columns of £ N. 




i 1 3K T 


Since rank [ijf] must be n, this condition gives 


n = rank 


30 , 


S N. 


i 1 3K T 


+ (n - p ) 
Lmax 


max 


= rank 


K 


o 


d<t> i 

i 1 <?K X t K 


(52) 


\ i 


Now consider the special neighborhood (\ , 6 q3 0), i. e., pole placement in 
the vicinity of the open-loop system. For this neighborhood, all partials 


aK 


0. (K) vanish save those for the first-order terms N^ . 

-*- 2 . jtj 


This gives 
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max 


K = 0 
o 


= rank j N J N* . . . N 1 . . . N r [ 
(12 ml m ) 


( 53 ) 


or verbally, the maximum number of arbitrary poles in the vicinity of the 
open-loop system is equal to the rank of the matrix of coefficient vectors 
formed from the system's numerators. If we assume that the system is 
completely controllable from at least one of its inputs, say u^, and that the 
matrix M(Q, y) has rank m, then the following bounds are readily established 
for p r 


max 


m < p 


max 


^ min (n, mr) 


(54) 


lK 0 = 0. 


Application to Systems with Compensators 


These results apply verbatim to systems which include dynamic compensators. 

til 

For example, suppose that q -order dynamic compensation is permitted in 
the controller L(z). We then append q integrators to the original dynamic 
system, i. e. 


x = F x + G^u + G 2 q 

x . 1 = x + u , .. 

# n+1 n+2 r+1 

X . 0 = X _lQ + U, c 

n+2 n+3 r+2 


x . = u , 

n+q r+q 


M(a y) x 


x 


n+1 


z = 


x 


n+2 


(37a) 


(40a) 


x , 
n+q 


62 - 



u 


k u kT 12 


K 21 K 22 


z 


(41a) 


Each integrator adds one input and one output to the original collection of 
inputs and outputs and thus adds additional N.(s) 0. (K) terms to equation (46). 
The new value p is again given by equation (52) and it follows from this 

ZXIdA 

equation that 


P 


max 


No 

compensation 




max 


th , 
q order 

compensation 


(55) 


This inequality assumes that the right-hand side is evaluated at 


K 11 K 12 
o o 


K = 
o 


K 21 K 22 
o o 


and the left-hand side at K = K 11 . If in particular we let K 11 = 0 and 
12 o ° o 

choose K q such the (m+l)^ measurement is connected to the control u 

(for which the original system is completely controllable), then an analogy 

to (54) yields 


m+q 


^ P 


max 


£ min [(m+q) (r+q), n+q] 

q L 1 order 
compensation 



0 K 


12 -i 


0 0 


(56) 
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Inequalities (54) and (55) hold in general for different orders of compensation. 
That is, p for q^^-order compensation equals or exceeds p for qj"* 1 - 
order compensation as long as q 2 £ q^. This fact can be used to determine 
the minimum compensator order required to place any prespecified number of 
poles, p, and in particular, the .order required to place all system poles 
(including the compensator). We simply start with q=l, check the rank con- 
dition (52), then go to q = 2, then q = 3, etc., until the rank equals n+q. 

The procedure will generate not only the compensator order but also a com- 
pensator design, via the function K(\) and 6 (\) -implied by equation (48). 


All this is best illustrated with an example: 


Consider the following system 


x = 


0 10 
0 0 1 

- d l - d 2 - d 3 


x + 


0 

0 

1 


u 


note: ji(l), . . ., i(r) | = ji(l) j = |3 | 


z = 


10 0 


0 1 O’ 


Suppose we want to assess pole placement capability without compensation 
and to design a minimum order compensator (if needed) for arbitrary place- 
ment of all system poles. 


Without compensation, we can use equation (44) directly to get the following 
first-order polynomial (in K) for the characteristic equation 
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3 2 

(S + PgS + PgS + Pj) = 

~s -1 0 "I f~s -l o 

(s^ + d^s 2 + ^2 S + + ^11 det 0 s -1 + K 12 det Os -1 

-1 0 0 0 -1 0 

or, in coefficient- vector -form 



Applying (52) we find that 



which means that only two poles can be assigned arbitrarily. The third 
remains beyond our control. Equations for all three can be obtained from 
equation (48) and the implicit function theorem, i. e. 


2 

specified roots: s + \ 2 s + 

unspecified root: s + 6-^ 
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Now let's explore arbitrary placement of all system poles via compensation. 
We begin by appending a single integrator, working up to two, three, or more 
only if required. The augmented system is 



with ji(l), ..., i(r)[ = ji(l), i(2)J = |3, 4 j 

Equation (44) now becomes a -second-order polynomial in K with first-order 
terms for multiplier functions 0^(K) = K l2 , K 21 , K 22 , K 23 and 

second-order terms for (K^K^ - K 12 K 21 >, (K n K 23 - K l3 K 21 >, and ‘(K^K^ 

K 13 K 22^' 

The result 
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with all other terms equal to zero. Note, in particular, that all first order 
coupling terms between system and compensator vanish, i. e. , terms for 
K 13 , K 21 , K 22 . Again applying (52)- gives 



CO 

eq 

W 

!___ 

0 

- K 21 

- K 13 

0 

- d l + K n " 

-1 

K 23 

" K 22 

0 

- K 13 

“ d 2 +K 12 

0 

-1 

0 

0 

0 

~ d 3 

0 

0 

0 

0 

0 

-1 


whenever f 0 (this can be shown by computing a determinant for columns 

2, 4, 5 and 6). Thus, arbitrary placement of all poles is possible with a 

first-order compensator in all neighborhoods except those about K Q such 

that = 0. This merely says that the compensator will do no good unless 

its output is fed to the input of the original system. [This corresponds to 
12 

the choice of K made to get equation (56). ] 


The compensator design and direct controller gains can now be obtained by 
computing the functions K(\), 5'(\), whose existence is guaranteed by the 
implicit function theorem. Since all poles can be placed, the function 6(\) 
is not present and the remaining function can be obtained directly from the 
polynomial P(s) above. This gives 


4 3 2 

specified roots: \{s) = P(s) = s + + x^s + x 2 s + Xi 

unspecified roots: 


K 23 

= 


K 12 

= 

d g - X 3 - K 23 d 3 

K 13 

i 

0 otherwise arbitrary 

K 22 

arbitrary 

K 11 

S3 

d l " ^2 ^23 d 2 + ^12^23 ~ ^ 

K 21 

= 

<-M - K 23 d l +K U K 23 )/K l3 


The final compensated system has the structure shown in Figure 14. More 
elaborate computational examples are treated in Appendix G. 
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Figure 14. Example Compensator Design 
-A Pole -Placement Algorithm 

In the example above, the system was sufficiently small to permit direct 
manual solution of the pole -placement problem. This will not be the case, 
of course, for large flexible -vehicle design problems.' To establish com- 
putational feasibility for the proposed quality measures, therefore, it is 
necessary to develop computerized solutions able to handle high order 
systems. These are provided by the Newton Raphson algorithm shown in 
Figure 15. 

This algorithm was used to carry out the .trial computations in Appendix G, 
which deals with 6th- and 17th-order dynamic systems. Computation .times, 
were very reasonable for both problems. For the 17th-order case with .two 
control inputs and four measurements, a single run consumes .approximately 
5 minutes on an. SDS 9.300 machine. About half of this time is devoted to the 
computation of generalized numerator coefficient vectors which present 
certain numerical challenges because of their large magnitudes -for high- 
order problems. At present, a generalized eigenvalue routine is used to 
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Figure 15. Pole -Placement Algorithm 
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perform the computations. This routine is discussed briefly in Appendix H, 
together with, an alternate procedure which may prove more efficient in 
•future applications of the pole-placement algorithm as an inner loop of the 
basic sensor .choice algorithm of Figure 13. The alternate procedure utilizes 
the fact that higher-order numerators of equation (44) can be expressed in 
terms of first-order numerators (Ref. 18). This fact significantly reducers 
the computations required fo evaluate the coefficient vectors for many 
sensor locations, y. 
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SECTION V 
CONCLUSIONS 


This study has produced partial answers to three fundamental problems in 
optimal controller design: (1) controller simplification, (2) sensitivity to 
model inaccuracy, and (3) sensor complement choice. Progress achieved 
on these problems are briefly summarized here along with unanswered 
questions. 

In the area of controller simplification, a gradient iteration method and' a 
compatible method for initialization were developed. With these methods a 
simplified controller was successfully designed for a rigid-body model of a 
launch vehicle. Parametric techniques were revised and used to design a 
simplified controller for a large flexible aircraft at one flight condition. 

The remai n ing question with respect to controller simplification is whether 
the techniques can be combined to handle flexible launch vehicle problems. 
There is no basic inconsistency between the methods, so it appears the 
question can be answered in the affirmative. However, verification can only 
come by actually exercising the techniques on a realistic problem. 

The major result with regard to model inaccuracy is that a boundary ‘point 
of the model's admissible parameters should be used to design an optimal 
insensitive controller. A necessary condition which such a boundary point 
must satisfy was derived. Computational aspects were considered and a 
numerical example corresponding to a flexible aircraft with unknown flexure 
frequencies was treated. The significant question remaining concerns 
admissible ranges of parameters. The magnitude of admissible parameter 
variations in the example was approximately 15 percent. At present such a 
figure must be determined experimentally on the computer for each individual 
problem. An a priori estimate is. highly desirable. 
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To gain .knowledge .concerning the quality of sensor complements, the .pole- 
placement capacity .of a controller with a given sensor complement was 
studied. The .following theoretical result was achieved for .multi -input 
systems: to determine the maximum pole -placement capacity .offered .by .a 
sensor complement, each input must be utilized, and indeed, the capacity 
of a given sensor complement is generally enhanced as mor-e inputs are • 
admitted. Furthermore, -this -capacity can be evaluated with .reasonable 
computational loads for compensated as - well as uncompensated systems. In 
fact, the minimal compensator required .to yield complete pole -placement 
capability can be determined computationally. Algorithms were developed 
and exercised on up to 17th-order examples. The question of how pole- 
. placement quality measures for a. sensor complement relate’to measures 
of quality with respect to controller performance has been left unanswered.' 

The positive results .-.aclneved in-each of the three -areas, of investigation make 
it possible to apply the stochastic constrained-response technology, to. the 
design, of a -controller -for mated .ascent of the Space Shuttle Vehicle. 
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APPENDIX A 

GRADIENT SEARCH FORMULATION 


The following flow chart describes the gradient method used for optimizing 
fixed-form controllers for time -varying systems. In the chart, K k (n) are 
the feedback gains of the k th iteration for time n, H k is the Hamiltonian of 
the k iteration, P (n) is the co -state matrix of the k th iteration for time n 

“■ ,v . > 

and X k^ is the covariance matrix of the k xn iteration for time n. 



A1 









If J = J** (the quadratic cost), the Hamiltonian is 
H k = TR [Q(N)H 1 (N)X k (N)H 1 ' :C, (N)] 

N-l 

+ TR C E At Q(n) (H 1 (nJ+D, (n)Kjn)M(n) )X (n-)’ (H^nJ+D^njK^n) 
n=o 

M(n)) T ] 


N-l 

+ TRC-S P. (n-+l)(X(n+l)-X(n))] 
n=o 

where Q(n) is the quadratic weighting matrix for time n, A-t is the sampling 
interval, M(n) is the measurement matrix for time n, and 1 H^(n), D^(n) are 
the usual system parameters. 

To solve for X k (n), use the following, sample- data covariance- solution' 
(forward integration): 

X k (n+1‘) = CA(n)+B 1 K k (n)M(n)],X k (n) [A(n).+B’ 1 K k (n)M(n)] T 
+(At)" 1 B 3 (n)W 1 (n)B 3 T (n) 

where. A(n), B^, Bg (n) , and- W^(n) are the usual system parameters. X k (0)< 
is known and is- constant. 

To solve for the. co-state matrices 1 P k (n), set 

SH k p k ( n+1 > - p k (n) 
a X k (n) At: 


A'2 



which has the solution (backward integration) with 
P k (N) = H 1 (N) T Q(N)H 1 (N): 

P k (n) = [A(n)+B 1 K k (n)M(n)] T P k (n+l)CA(n)+B 1 K k {n)M(n)] + 

At [ H 1 (n) +D 1 (n) K k <n) M(h) 3 ' T Q (n) [ H 1 (n) +D ± (n) K k (n) M(n) ] 
These equations yield the solution 

c_. 

* H k 

7 = 2B 1 T P k (n+l)[A(n)+B i (n)K k (n)M(n)]X k (n)M T (n) 

K k (n) | 

+2 At D 1 T (n)Q(n) Ch, (n)+D 1 (n)K k (n)M(n)]X k {n)M T (n) 

with 

AH, 

k 

= 0 . 

SK k (N) 


To choose AK k (n) small such that 


N aH k 

AJ = E Tr [ AK, (n)3 


n=0 


SK k (n) 


0 , 


we choose 


3H, 


Tr [■ 


AK k (n) 


AK k (n)3 £ 0 
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for all n, which* is true if 


3H, 


9K, fa . 4K k < h > * 0 


3H, 


for each i, j, and n, which is true if 


oK 


k. 


ii n) 


and AK, (n)' are of opposite sign, 
ij 


This was accomplished by setting- (for e > 0) 


dH k 

AK, (n) - - e|K (n) | 

ij xj SK k (n) 

ij 


£ / 3H k 

U j,n 


3K, (ri), 

k i;1 


If it is desired that J'* (the upper bound on the probability of mission failure), 
is minimized instead of J**, simply set all the Q(n) of the above equations to- 

3J* 

Q(n) = Q (n) = 

R as k (ir) 

where 

S k (n) =■ [H 1 (n)+D 1 (n)K k (n)M(n)]X k (n)[H 1 (n)+D 1 (n)K k (h),M(n)] T 


t^or the launch vehicle, we have a single input; this means i=l. Normalizing 
over each gain individually consists- in removal of summation over j in' 
radicand. 


A4 



This is true because 


= Tr [• 


5 Jr 


3S k (n) 


3K, (n) 


3S k (n) 


3K k.. (n) 


-3 




where 


3 j** 


£>K (n) 

k ij. 


3S k (n) 

= Tr CQ(n) ] 

&K k (n) 

ij 


which are identical if 


Q(n) = Q, (n) = 




^S k (n) 
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APPENDIX B 

SIMPLIFICATION OF TIME -VARYING GAINS 


A method is proposed for further simplification of time-varying optimal con- 
trollers through time linearization. It is proposed that the gain matrix K{t) 
be constrained to be piecewise' linear with respect to time. That is, let 

K(t) = a • (t-t V + b- for t , s ts t , & = 1, P (Bl) 

where t = 0 and t = T. 
o p 

The method considers optimizing via a gradient scheme with respect to the 

a , and b . with a constraint for continuity of K(t).- The breakpoints- 1 , are 

& A ... . , 

determined visually from, practicalized' gains obtained prior to these param- 
eter optimizations. 

Also discussed is. the determination' of the deterministic input for the 
simplified controller'., 

FORMULATION; OF THE PROBLEM. 

Figures Bl through-B5- show how five- time-varying gains may be split up in a 
piecewise time.- sense.. If seems’- that this could be done- reasonably with six 
breakpoints;, t through f g .- Actually, the- set of breakpoints would be’ the union 
of the- necessary breakpoint’s' necessary for each individual time -varying gain. 
The knowledge- of the: physics^ involved along- with the’ visual observation of the 
gains would: determine:- where' the: breakpoints- should be. From a practical 
point: of view it, Is' advantageous: to- have as’ few breakpoints as possible. 


Bl 




6100 8-00 tfl-OO 

TIME SEC CX10 *J 


Figure Bl. Optimal $ Gain with Initial Approximation 
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6.00 8.00 10.00 
TIME SEC (X10 *1 


Figure B2. Optimal z Gain with Initial Approximation 







6.00 8.00 10.00 
TIME SEC CX10 1 ) 


Figure B3. Optimal j3 Gain with Initial Approximation 
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Figure B4. Optimal 0 Gain with Initial Approximation 
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Figure B5. Optimal z Gain with Initial Approximation 


B6 


Given the set of breakpoints, we then define an initial set of gains (in dif- 
ference equation formulation) 

K Q (n) = a ^ (n-n^) +b^ ; <; n <. l-l, p ' (B2) 

where 


b = K*(N)t (P.racticalized final time gain) 

Po 

\ = Vi 0 ( VVi ,+ Vv 1=1 H 


(B3) 

(B4) 


and the a are the slopes of the lines determined graphically or through 

■Vq 

parameter optimizationf where in conjunction with equations (B3) and (B4) 


o 


I 

n=n 


r 


i-i 


b „ - K*(n) 

l o 



(n-n )“ 
& 


(B5) 


l-l 


tb naay also be determined through parameter optimization by solving the 
Pq following two equations simultaneously: 


o 


N 

2 

n=n 


p-1 



r -1 

) 


b - K*(n) 

(n-N) ( 


Pq 


1 


) 


N 

£ 

n=n. 


p-1 L 


N 


K*(n) - a (n-N) 
P r > 



(n-N)' 


/( N -n p _i) 


,2 .. 


in which case j [K*(n) - K (n)] is minimized with respect to a 

n=n P -i p ° 

and b 

p o 


B7 



n 4 

where K*(n) are the set of practilized time-varying gains. Then Z 

n=n 

2 

[K*(n) - K (n)] is minimized with respect to a . Time n is equal to N; 

° ■ z o p 

n Q is equal to zero. 


4-1 


When optimizing on the a^ and with the constraint for continuity of K(t), 
the Hamiltonian is 

H k = TR[Q(N)H 1 (N)X k {N)H 1 T (N)] 

fN-1 
S 

n=o 


+ TR <j Z AtQ(n)[H r (n) + D 1 (n)K k (n)M(n)]X k (n)[H 1 (n) 

( n=o 

t) ' f N_1 ' ) 

+ D 1 (n)K k (n)M(n)] 1 |+ TR j Z P k (n+1) [X(n+1) - X(n)]| 


(B6) 


where Q(n) is the quadratic weighting matrix for time n, At is the sampling 
interval, P. (n) is the costate matrix and X T (n) is the cpvariance matrix of 
the k n iteration for time n, M(n) is .the measurement matrix for time n, and 
Hj(n) and D^(n) are the usual system parameters. The K k (n) are defined as 


K k (n) = a^ (n-n^,) + b ^ n ( j_ 1 ^ n 


(B7) 


where 


b = K*(N)t 
Pk 

with the constraint that 


(B8) 


b 4 k = a ^+l k <n 4” n A+l > +b 4+l k ; A 1 '*** jP 1 


(B9) 


fSee footnote on page Bll. 



which can be written as 


b 




P 

S 

j=j&+l 





+ b 

Pk 


(BIO) 


From equation (BIO), we see that the b , are functions of the 

A k 

is greater than t , and b , and not the a . 

Pv j8i. 


hi 


9 V* 


To solve for X^(n), use the following sample data covariance solution: 


X k (n+1) = [A(n) + B 1 K k (n)M(n)] X^n) [A(n) + B 1 K k (n)M(n)] T 
+ (At)" 1 B s (n) W 2 (n) B 3 T (n) 


(Bll) 


where 


A(n), Bg(n), and W (n) are the usual system parameters. 
X k (0) is known and is constant. 


To solve for the costate matrices P k (n), set 
P k (n+1) - P k (n) 

^X k (n) At < B12 > 

which has the solution with P (N) = H (N) T Q(N) H (N) 

k l 1 

p k (n) = [A(n) + B K (n)M(n)] T P (n+1) TA(n) + B K, (n)M(n)] 

K 1 K (B13) 

+ At [H^(n) + D^(n) K k (n)M(n)] T Q(n) [H^n) +D 1 (n)K k (n)M(n) ] 


B9 



These equations yield the solution 


dH^ N 


da « k = n=o SK k <n> 


dK k (n) 
da . 


dK{n)' 


d a 


0 

, n — n . 

& 

11. 1 ~' n A 

£,~L & 


l -l 
n s n 


sn 

^•n. s n„. 


a 


j6-i 


(B14)' 


where 

5H k 

9K k (n) 


2B 1 T P k (n+1) (TA(n)+B- 1 .(n) K^nJMfn) ]' X k (ri) M T (n) 

rp 

+2 At D 1 T (n) Q(n) QHjtiu) H-D^n) K k (n)M(n)] X k (n) M (n) 


(Bl-5). 


Note that 

dK(n) db.. 

= - 3 -J- == n- -n., for' n'. . £>n. <. n., j < 4 
d a^ da^ jJ-1 4- 3. - ' 1 ' 3’ 

We can, use equation (BIS’)' computed backward in time to choose AK k (n) in 
a gradient search,- where. 


AK k (n) - Aa^ (n-n^) + Ab^. y n^_ I £ n s l- l, 
’ k k 

where. 


(-B16-) 


A a., = -e- 

4 kr 




ki‘ 




a. 3 .-. 
4=^ P- V *k 


T 5 a-- 


\i 


(B1‘7) 


th* 

The subscript i represents, the- i component' and e is- the maximum 
percentage change- in the a: .. 


BIO 



Also 


Ab 


A k-1 

P 

= z Aa. (n. 
j-A+l J k 


-n j) + < b p 


k 



where 

b = K*(N). t 
P k 


(B18) 


The gradient search is set up as in Appendix A. If it is desired to minimize 
a nonquadratic cost <T : , let 


Q(n) = Q k (n) = 


dS k (n) 


(B19) 


+ This b 

P k 


may also be a variable. 


where 



N 

nlc ^ 5 


dK k (n) N a H^ 
®p k = n=o 


instead of being fixed since this parameter is not constrained. 
Note 

dK, (n) 

= 1 for all n from (B7) and (BIO) 

Pk 


Bll 



where 


S k (n) = [H x (n)+D (n) K k (n)M(n)] X k (n) [H^nJ+D^n) K k (n)M(n)] T 
as discussed in Appendix A. 


THE DETERMINISTIC INPUT 

From the perfect .sensing gains K (n), the deterministic input f (n) is 
defined. The deterministic controller u is then 

u(n) = K (n)x{n) + f n (B20) 

Since this is the optimum deterministic controller, it may as well be used 
with the simplified controller gains K(n), and a new deterministic 

a- 

input f ’ (n), where 

f^ (n) = u(n) - K(n)M(n)'x(n) (B21) 

where u(n) is the controller defined in equation (B20) and x(n) is the mean 
state vector found by 

x(n+l) - A(n)x(n) + B .u(n) + B„(u)v (n) (B22) 

1 & cu 

which is the same as the ; optimum since u(n) is optimum. 

In order to further s.implifydhe deterministic input to the form 

f(n) = c * .(n-n •) + d for n. - ^-n £ n, ; Z~l, . . ., q 
Z Z Z x>~ l 


(B23) 



an initial set of c and d can be found with a straight parameter optimiza- 

i a 

tion where 



f (n) - c (n-N) 
L q ° 




f*(n)3 


(n-n .) 
I 


< N - n q-l ) 



(B24) 


, )2 

2 (n-n ) ; 1=1, 
n ~ n z-i 


-, q 


and 


\ “ Vi Q (n i “i+i* + d ^ +1 0 J x=i, 


(B25) 


and then use a gradient method to improve on these. It is obvious from 
Figure B6, that the same breakpoints as for the K(n) cannot be used. For 
the gradient method we would use the same procedure as for the K(n) 
except that the gradient is now for the k iteration 


a H . * 


JlIL 


3 C 


K n=n „- 


je-i 


aMn) 


N 3 h ay*) 


S (n) ac ( 
n=o ° k ) 


where 


Sf k (n) 


o 

n-n . 

& 

n . -n . 

j &-1 l 


n^n 

Vl S “ n 4 / 


n s n 


a - 1 
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Figure B6. Optimal Deterministic Input 
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and 



N 

E 

n=o 


5H 

3 f k (n) 


df k (n) 



N 

2 

n=o 



afjj.(n) 


since 


Sf k (n) 



1 for all n. 


where H is the Hamiltonian associated with the mean responses and 
dH B T 

= -Jj- x k (n+l) + 2D 1 T (n)Q(n) jtH^n) + D^n) K(n)M(n)] \(n) 
k 

+ D 1 (n) f k (n) + D 2 (n) v^(n)[ 

where, 

\ k (N) = 2H 1 (N) T Q(N) [H 2 (N) x k (N) + D 2 (N) v^N)] 

f rn 

X (n) = [A(n) + B : K(n)M(n)] T x k (n+l) + 2At[H 1 (n) + D^n) K(n)M(n)] T Q(n) 
| [H 1 (n) + D^n) K(n)M(n)] x^n) + D^n) f fc (n) + D 2 (n) vjn) j 

and 

x k <n+l) = [A(n) + B 1 K(n)M(n)] 5^(n) + B 1 f k (n) + B 2 (n) vjn) 
with x^O) = x^. 

The procedure is then to compute x k (n) forward, and x k (n) and the gradient 

backwards, compute c and d by 

^k+l q k+l 


B1'5 



lc ,i :C d 


,d =*4_ + ;Ad 

q .k+l q k % 


.or compute W (I,,> ;b - y 

■ f k+-l^ n ) = - f k ( ' n) + - A ^k^ 


■Ayn) --= J& - >(n-n^) -+ ,A<y ; ; ; n^ ^ ^ n.<: , £=% <1 

;k * k 

The new .are used ,to repeat tt he -procedure for the .next step. 
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APPENDIX C 

TIME -VARYING GAINS FROM GRADIENT ITERATIONS 

This appendix consists of graphs of the gains found in Section II, Iteration 
Method for Time-Varying Systems. They are presented as the following 
figures: 

Cl. 0 Gain, Iterations from Perfect-Measurement Gain 

02. z Gain, Iterations from Perfect -Measurement Gain 

C3. p Gain, Iterations from Perfect -Measurement Gain 

C4. 0 Gain, Iterations from Perfect-Measurement Gain 

t 

C5. z Gain, Iterations from Perfect -Measurement Gain 

C6. 0 Gain, Iterations from DC Approximation 

C7. z Gain, Iterations from DC Approximation 

C8. p Gain, Iterations from DC Approximation 

C9. 0 Gain, Iterations from DC Approximation 

CIO. z Gain, Iterations from DC Approximation 

Cll. *0 Gain, Iterations from Averaged Approximation 

C12. z Gain, Iterations from Averaged Approximation 

C 13. p Gain, Iterations from Averaged Approximation 

C14. 0 Gain, Iterations from Averaged Approximation 

C15. z Gain, Iterations from Averaged Approximation 

C16. 0 Gain, Initial and First Iterations 

Cl 7. z Gain, Initial and First Iterations 

Cl 8. p Gain, Initial and First Iterations 

C19. 0 Gain, Initial and First Iterations 

C20. z Gain, Initial and First Iterations 
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Figure Cl. 0 Gain, Iterations from Perfect-Measurement Gain 
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Figure C2. z Gain, Iterations from Perfect-Measurement Gain 
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TIME SEC CX10 l ) 


Figure C3. /3 Gain, Iterations from Perfect-Measurement Gain 
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TIME SEC CX10 

Figure C4. 0 Gain, Iterations from Perfect-Measurement Gain 
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Figure C5. z Gain, Iterations from Perfect-Measurement Gain 
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Figure C6. 0 Gain, Iterations from DC Approximation 
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Figure C7. z Gain, Iterations from DC Approximation 






















6. DO 8.00 10.00 

TIME SEC CX10 l l 


Figure C8. (3 Gain, Iterations from DC Approximation 
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6.00 8.00 10.00 12.00 11.00 16.00 
TIME SEC (X10 


Figure C9. 0 Gain, Iterations from DC Approximation 
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Figure Cl 1. 0 Gain, Iterations from Averaged Approximation 
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Figure 02. z Gain, Iterations from Averaged Approximation 
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Figure C13. |3 Gain, Iterations from Averaged Approximation 
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Figure C14. 0 Gain, Iterations from A 
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Figure Cl 5. z Gain, Iterations from Averaged Approximation 
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TIME SEC 1X10 


Figure Cl 6. 0 Gain, Initial and First Iterations 
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Figure Cl 8. /3 Gain, Initial and First Iterations 
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Figure Cl 9. 0 Gain, Initial and First Iterations' 
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Figure C20. z Gain, Initial and First Iterations 
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APPENDIX D 

COMPUTATION OF SENSITIVITY COEFFICIENTS 


In this appendix three computational methods of determining sensitivity co- 
efficients of the performance index for parameter variations of the system 
are presented. 

Let the system equation be 

x = F(p) x + G(p) u ; x<o) = x q <D1) 

where F is a nxn matrix, 

G is a nxr matrix, 
p is an m-vector of parameters 

Let the performance index be 

CO 

J = J (x T Qx + u T Ru)dt (D2) 

o 

where Q is a nxn positive definite matrix, and ft is a rxr positive definite 
matrix. Let the control law be 

u = Kx <D3) 

where K is a rxn feedback gain matrix. Then the optimal feedback gain 
matrix K* which minimizes J, becomes 

K* = -R' 1 G T (p)P < d4 > 

D1 



where P satisfies the following Riccati matrix algebraic equation 

-PF - F T P + PGR _1 G T P - Q = 0 (D5) 

and the optimal performance index J* may be written as 

. J* = { x o T Px o } = Tr jp X q | (D6) 

T 

with X denoting x x 
o & o o 


On the other hand let us assume that K is simply a feedback gain matrix which 
stabilizes the system <D1). Then the corresponding system equation becomes 

x = (F + GK) x ; x(o) = x q (D7) 


and the performance index is 

,= L T fe <F+GK) T t (Q+K T RK)e (F - t<3K > t ]x„ dt 
JO o 

o 

= x o T [J e (MK ) Tt (Q +K T RK)e (F « 3K>t dt]x 0 
o 

= x T Ax = Tr )a X } 
o o ,< o> 


(D8) 


(D9) 


where 

A & J e <I--K3K) T t( Q+K T RK ) e (F4GK)t d t 
o 


(DIO) 


and it satisfies 


(F+GK) T A + A(F+GK) + (Q+K T RK) = 0 


(Dll) 


D2 



Note that P is a function of only the parameter p, whereas A is a function 
of both the parameter p and K. 

Now we can define the performance index error due to using a simple stabili- 
zation matrix K rather than the optimal feedback gain K* for any value of the 
parameter p as follows, 

e(K,p) = J(K, p) - J*(p) 

T T 

= x Ax - x Px 
o o o o 

= Tr | (A - P) X j 

- Tr[ |a(K, p) - P(p)J X Q ] (D12) 

-1 T 

Note that if K - K* = -R G P, then A equals P and e(K, p) is zero. 

The function e(K, p) can be expanded into a Taylor's series in terms of K and 

p at any nominal values of K and p as: 

0*0 

e(K,p) = e(K 0 ,p o )+V K e[K-K 0 ]+Ve[p-p o ] 

+ \ [K ' K o ]T V Kp e[ P^‘o ] + i [ P-Po- Ty pK e[K ~ K o ] 

+ itK-K o ] T ^[K-K^+iCp-p/ V p e[p-p o ] 

+ higher order terms (D13) 

where K and [K-K^lare considered as vectors and all partial derivatives are 
evaluated at (K Q , p Q ), V^e is the- (rn) row vector 



D3 



and V T ^_ is' the (rn)'x(m) matrix 


e = 
Kp- 


2 

3 e 


2 

3 -.e 


2 

3 e 


-2 
3 -e 


2 

3 e 


9K H 5P i 3' K li3P 2 ”* SK 11 SP m 


"2 
3- e 


3^e- 


SKg j. 3P^ 3 Kg y 3P2 3 K 2 -^3P 


m 


2 

3 e. 


2 - 
3 e 


& K rn aPi a K rn ap 2 :'' aK rn ap. 


m 


Similarly other second-order partial matrices are defined. 


Since e =■ JfK, p)' - J*(p) 

- Tr j A(K, -p) - P(p) } X 6 ], . 

_ 3 £ _ aJ _ 3 J * - 

3k‘ 3k 3-k 

=■ Tr K|# - HX f 

= Tr X Q | where, k denotes an element of K (D14)'’ 

Thus, in order to obtain a Taylor's series- expension of e, we must obtain 
the partials of A. and Pwith respect to K and p. These partials are called 
. sensitivity, coefficients- of performance. indices. 

Since.- A and P are. implicit functions of K and p, there is no direct method 
to determine 'these, partials. The 'partials can be determined by solving 


D4 



Lyapunov type matrix equations generated by applying the implicit function 
theorem to (D5) and Dll). Sample sensitivity coefficient matrix equations 
follow: 


(F-H3K) T (^-)+(^-)(F-K3K) = -jCAG+K T R]{E ij )| T 

- j [AG +K T R ] (E. ■)} (D 1 5) 

where E. . denotes an (r)x(n) matrix with the ij** 1 element equal to one and 
all other elements zero. 

( f+gk) t (-|^) +(-^-)<f+gk) = -(|E- + ||-k) T a 

' + lf K ) (D16) 

(F4GK*) T (^)+(^)<F-K3K*) = + P 

p/SF ] cG K*\ (D17) 

} 

Similarly, second-order sensitivity coefficient equations are of the form 


D5 



rn/ O \ I O «*'-»■ l 

(F+GK) " " L. / + Uf Tk j (F+GK) 

' ^ij^km' VaK ij aK km / 


#^-)g(e ) 

jK krr/ 1J . 


■ 3 rr) G «W 


BA 


,BK 


ki 


)g (E ij> 


13 


"IT 


(E..) R(E, ) 
13' km' 


3 K. . j G ^ E km ) 
^*3 


< E l;j> R < E km> 


(F+GK) 




= - (GE..) T GE-. 

i] Sp k Bp k 13 


3F , ^G tA BA _ cA / 3F 4 . 3 G xr 
1 BK,, " 5 K.,UP k SP k . 


, 3 P V BP 


k 


13. ■ 13 


-^-,e...\ ■ A -a|^-e;.., 
BP k 13/ . BP k 13; 


Three: methods- for" computing, these partials are now discussed. 


ADJOINT’ METHOD-’ 


Let the adjoint matrix-' equation, he 


(F'+GK)'W + W-(F+GK) =- X r 


(D18) 


(D i 9 )' 


(D 20 ) 



Since (F+GK) is a stable system matrix and X Q is a positive semidefinite 
matrix, W is a negative semidefinite matrix. 


Consider 




aKij 


Fy °! 


(D21) 


Substituting (D20) into (D21) produces 

^§77 = Tr (^— )C<F+GK)W + W(F+GK) T 3 




(D22) 


Using the trace properties, Tr(AB) = Tr(BA), Tr(A L ) - Tr(A), we obtain 
-|5— = Tr |[(F-K3K)W + W(F+GK) T ](-^— )| 

j[<F4CK>w(JC) + w(F4aK) T (^-) | 


= Tr 


= Tr (F +GK)W + (F-HBK) T '(-^-)w] j 

( ij d ij 

= Tr |[( F 4GK) T (-^-) + (^§ 7 ) (F+GK) ] W ' 


The expression inside the bracket of (D23) is the left side of (D15); hence 
c>e 




- -Tr 


|[(AG+K T R)(E i;j ) ] T + [ (AG+K T R)(E i;j )]| W 


(D24) 


As we can see, the right side of (D24) is known and hence - ■ can be deter - 

3K ij 

mined for all i, j; i = 1 , . . . , r; j = 1 , . . . , n. 


D7 



Define another adjoint matrix. V by the equation 
(■F+GK*)V + V(F+GK*) T = X Q 

Then 


Tr 


~SP 

' — • X 

BP k °' 


= Tr 


~gP gP 

— (F+GK*)V + V(F+GK*) 

_3P k ‘ 5P k 



c)P' . - .T a p ■ 

— - (F+GK*)V + (F+GK*), — V ■ 

sp k ap k 


= Tr 


- .c)P v q 1 $P 

!■ (F+GK*) + (F+GK*) 

1 '*P k . . - . 5P k . 

. — 



(D25) 


(D26) 


Using. (D-17;) in. equation (D26) yields- 


Tr 


a-P- 


- X 


.Wk-. 


o' 


= -Tr 


• /5 f- 5 g v t /SF b g 
— +■ . K* jV P + P" [ + 

v¥k. •• sp k J 


K* 


^P k ^ p k 


Vl 


Therefore 


ae- 

/d A 

a p V .. 1 

= Tr 

f— " 

)' x 0 

3P k _ ^ 

:\^ k 

a? J. . 


= -Tr 1 


a-F aG; v 
— - + — . K 


■ T' 


/aF- s g \ 
A + A + K' 


\ 5 p k - 


W 


+Tr' 


^ SPfc ?>P k 

79* ?>G; V /SF S G 

I + K-J, P + P +' Kp 

W a p k . r W k ap k 


T 


• VT (4027)' 


D'8 



The right-hand sides of (D24) and (D27) involve only known terms and V and 
W. Thus the first-order partials can be computed using (D20), (D24), (D25) 
and (D27). From equations (D24), (D26) and (D27) the following properties 
can be established: 

Property 1 : e(K, p ) = 0 

K=K*(p o ) 

Proof: K-K* implies A = P since the equations (D5) and (Dll) are 

identical for K=K*. Therefore 

e(K,p o ) = Tr[(A-P)X o ] =0 

K=K*(p ) A=P 

1 o 

Property 2: Vj^ e = 0 

K=K*(p) 

Proof: From (D24), 



D9 




Equations for the second-order partial- derivative matrices may be obtained 
using the adjoint matrices and equations of the form’ of (D18) and (D19).- .For 
example 



DIO 



,A 


2 

> e 


As we can see from (D28) we need 


only way we can obtain 


2 

3 e 


c)K.- 


to determine 


aK..aK. 

° l] 0 km 


The 


3K-. 

il 


is to solve the matrix equation (D15). We 


v 

must solve rn matrix equations to have V T ^ T ^e. This is a severe drawback 
for the adjoint method of determining the partials. 

From (D28) we can obtain a direct proof of the following property of the 
second partial. 


Property 4 : V KK; e 


is a positive semidefinite matrix. It is 


K=K* 


positive definite if X Q is positive definite. 


Proof: 




Thus- 


5 k 


If K=K*, then the right-hand side of equation (D15) is zero. 
= 0 where k denotes any element of K. Hence, from <D28) 


K=K* 


2 

a e 


° ij° km 


' - 2 Tr[(E..) T R E km W] 


K=K* 


Direct calculation reveals that TrlE.. 1 R E. W] = r.. w .. The matrix 

L i] km J ik mj 

W is symmetric so that r„ w = r.. w. . Using the Kronecker product 
J ik mj ik jm b 

notation we can then write: 


V KK e 


= -2 RXW 


K=K* 


Dll 



r is assumed to be positive definite and W is negative semidefinite (definite 
if X Q is definite). The eigenvalues of the matrix RXW are X^ where x i 
are the eigenvalues of R and p. are the eigenvalues of W (Ref. Dl). Since 
R and W are symmetric, x^ > 0 anc * Mj s 0 for each i and j. Thus X^Mj ^ 0 so 
the eigenvalues of -RXW are greater than or equal to zero. Strict in- 
equalities follow from the assumption that X q is positive definite. 

Next we will consider the second method of computing the partials which is 
a vector equation method. 


VECTOR EQUATION METHOD 

The major drawback of the first method is that solutions of many Lyapunov- 
type equations for different inputs (right sides) are required. Since the left 
side of the equation does not change, we may not need to solve the Lyapunov 
equation for each different right side. There have been numerous investi- 
gations which converted the Lyapunov equation into linear algebraic equations 
(Ref. D2, D3). We will use Bingulac's technique (Ref. D3). 

Let the Lyapunov-type equation be 
PA + A T P = -Q 

Then Bingulac showed how to convert this equation into the following linear 
algebraic equation: 

BP = -Q 
v v 


D12 



where 


P v = ^ P n’ P 12' * * *’ P ln 3 P 22’ * * * J P 2n J P 33* * * ** P n-1 n-1' P nn? 

T \ 

\ = *^12’ ‘ * *' q ln’ q 22 5 ‘ • ’* q 2n' q 33’ *** ,q n-l n-T q nn' 

B = jb. . | is a mxn matrix depending on the elements of A and 

n(n+l) 

m . 


Let us apply this technique to (D15). 

(F+GK) T / — \ +(— ^ (F+GK) = -C 

where C = { [AG+K T R](E..)} T + { [AG+K T R]{E..) j. The corresponding 
linear algebraic equation is 


D 




(D29) 


n(n+l) 

where D is a 

2 

of (F+GK), and 


n(n+l) 

matrix and its elements depend on the elements 

, 2 




v 




9 A 12 

a A ln 

^22 

^ A 2n 

^33 

^A 

nn 

> • • • 9 

a K., 
d i] 

^K. . 
° 10 

S K ij‘ 

• • • 9 ) 

aK.. 
d 13 

j • • • J 

SK.. 
o* x;] 

S K. . 
d i] 


< = (C H’ C i2' * * *' C ln > C 22' C 23 J * * " C 2n' C 33' * * *' °3n J * * * J C n-1 n-F 


C ). 
nn 
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The equations corresponding to second partials of A may be expressed in 
the form of (D29) with the same D matrix. The right hand side (the C^'s) 
would be linear functions of the first partials. Thus the second partials could 

be expressed in the form 
D -1 [T 1 D" 1 T 2 + T g ] 

where the T. are well defined matrices and vectors of appropriate dimensions. 
The elements of T\ depend on F, G, K, Q, R and the indices of the variables 
with respect to which derivatives are being taken. Thus, the second partials 
can be easily computed once D _1 is determined. However, this method has 
a different drawback in comparison with the first method, namely unreasonable 
computer memory requirements. For example, let the order of the system be 
20. Then the number of elements of D is [20(20+l)/2] 2 = 44, 100. This 
storage requirement is too severe. 

The first method for computing the partials requires solutions of many 
Lyapunov equations, whereas the second method requires extensive computer 
memory. A third method, described below, is a compromised version of 
the two previous methods. 


POWER SERIES METHOD 

Jameson (Ref. D4) formulated a method of solving a Lyapunov-type equation,, 
using a power series of A up to n th power, without increasing the dimensions 
of the equations. 

Let a Lyapunov-type equation be 
AX + XA T = C 


(D30) 



Then X is expressed as 


X = Z 


-1 


r s ’ 1 ’ 1 <-i) i+j a. (A T ) < i ) 

Li=0 j=0 1 


where 


(D31) 


Z = 2 ( — 1 ) 


n 


r n/2 
Z a 
i=0 


n-2i 


. A 


2i 


; A° = I 


n 


n 


Xl - * 1 

Z a. TrCA 11 " 1 ) 
i=0 1 


Applying (D31) to 


cA 


cA 


, we obtain 


km 


cA 


= Z 




km 


Z ^ E X (-l)^a. (F+GK) (n i_1 ^ /-1(E ) T [AG+K T R] T 
i=0 j= 0 1 \ km L 

T (j) 


- [AG+K X R] (E km )| (F+GK) 


(D32) 


where 


Z = 2<-l) 


n 


n/ 2 

y a (F+GK) 
r A n-2i 
1=0 


2i 


; (F+GK)° = I 


a = - 
n 


1 

n 


n-1 

E a. 

i=0 


. Tr | (F+GK) 11 " 1 1 
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Note that Z" 1 and a are fixed quantities for a given (F+GK). Then 


n 


Be / . 

= Tr ( X 




km 




O 


m 


= Tr [Z 


-1 


r s" 1 ' r S 1 " 1 (-l) i+j a. (F+GK) (n “ i " 1 ^>/_(E km ) T [AG+K T R3 T 
i=0 j=0 1 


(j) 


-[AG+K T R](E km )| (F+GK) T ]X q 


(D33) 


Using the trace properties, we can simplify (D33) as 


3 e 


= -TR 


BK. 


km 


n-1 n-i-1 ... T (j) _i 

2 £ (-1) 3 a. ((F+GK) ) X q Z (F+GK) 

i=0 j=0 1 


(n-i-l-j) T 


^n-1 n-i-1 ... m (j) 


+ E E 
V i=0 j=0 


(_l) i+ 3 a . ((F+GK) T ) X o Z“ 1 (F+GK) (n " 1 ~ 1 3) 


(AG+K T R) (E km ) 


= -Tr[(Y T +Y) (AG+K T R) ( E km >3 (D34) 


where 


n-1 n-i-1 


(j) 




II- J. ll-I-I - , • rp - 1 

s s (-1) 1 3 a. ((F+GK) 1 ) X Z "(F+GK) 
■ i=0 j=0 


Note that (Y T +Y) (AG+K T R) is a fixed quantity for a given (F+GK). 


The second partial," 


2 

9 e 


B E ij B E km 


is expressed as 


D16 



2 

a e 


J A 


^ij^kzn 


= Tr 


X 


• SV^m 


o 


•Tr 


(Y T +Y) 


aA 


aK 


km 


|G (E..) 


-Tr 


(Y 1 +Y) 


aA 


aK, 


G (E. ) 

km 


i] 


-Tr 


(Y T +Y) (e..) t R (E. ) 

ij km 


Hence 


2 

a e 


aK aK. 

° rs° km 


(D35) 


= Tr 


T -1 
(Y +Y)Z l 


n-1 n-i-1 . . , . . 

E 2 (-1) 1 3 a.(F+GK) ln 1_ 3 TE. ) (AG+ 

U=0 j=0 1 km 


*<j> 


K T R) T (F4GK) T ' I G(E ) 

X k) 


j)" 


+ Tr 


T -1 
(Y i +Y)Z 1 


(F+GK) 


E Z X 1 ("l) (J) a. (F+GK) (n " i_;L j, (AG+K T R) (E, ) 

Li=0 j=0 1 km 


,(j) 


G(E ) 
rs 


+ Tr 


(Y T +Y) -1 
Z 


2 2 1 (-l) 1+;i a.(F+GK) (n “ 1_1 "^(E J T (AG+ 

L i=0 j=0 1 


rs 


rp rp rp(j) 

K R) 1 (F+GK) jG(E km ) 
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+Tr 


{Y T +Y)Z 1 


rn-1 n-i-1 


E T (“l) 1+ ^a. (F+GK)^** 1 " 1- ^ (AG'+ 

L 1=0 j-0 1 


K T R) (E rs )(F+GK) T 


(j)“ 


G(E km> 


-Tr [<Y T +Y)(E rs ) T B(E km )] 


(D36) 


As we can see from (D36), it is required to compute (F+GK)^ to determine 
the second partials. If we .store all powers of (F+GK), it requires 20 x 400 = 
8, 000 words memory for a 20th-order system. On the other hand, we can 
compute (F+GK)^ from (F+GK) stored in the memory which requires only 
400 words memory but this requires a large amount of computing time. 
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APPENDIX E 

PARAMETER VARIATIONS EXAMPLE 


f 

Numerical results obtained in the study of sensitivity to parameter variations 
for the C-5A example are presented in this appendix. These results demon- 
strate the applicability of design of optimal insensitive controllers to a real- 
istic problem with significant range of admissible parameters. The limita- 
tion in the magnitude of the parameter range is apparent for this problem. 

The model of the vehicle used in this study is described, followed by a de- 
scription of the problem treated. Numerical results indicating the nature of 
dependence of performance on parameter variations is then summarized. 
Finally an analysis of eigenvalue dependence on parameter variations is pre- 
sented which indicates limits on the admissible range of parameter variations 
for which the proposed approach to insensitivity is valid. 


MODEL DESCRIPTION 


The longitudinal axis of the C-5A, a large transport aircraft, was chosen for 
the constant coefficient example described in Section II. This vehicle was 
chosen because of the inability in a previous study to derive a simplified con- 
troller from an optimal controller. Also the data was readily available. The 
single flight condition treated is the low-altitude cruise condition with 50 per- 
cent fuel and 50 percent cargo indicated in Reference 10. The original data 
available consisted of three rigid-body modes, fifteen symmetric structural 
flexure modes described in Table El, lift-growth effects expressed in terms 
of Wagner and Kussner functions for the wing, horizontal tail, vertical tail 
and fuselage and gust penetration effects represented by time delays. A 
design model was derived for the optimization study reported in Reference 10. 
This model included two rigid-body modes {vertical velocity and pitch rate). 


El 



structural bending modes 1, 3, 4, 6, 7 and 11, Kussner lift-growth modes, 
a second-order wind filter and three first-order actuators. The expected 
value of a quadratic form in stresses and stress rates at two wing stations and 
three fuselage stations, including the horizontal tail root, and normal accel- 
erations at four fuselage stations was chosen as the performance index. 
Aileron, inboard elevator and spoiler control surfaces were used for control. 


Table El. C-5A Mode Descriptions 


Mode 

Frequency 

(Hz) 

Description 

1 

0. 750 

First wing vertical bending 

2 

1.780 

First wing chordwise bending 

3 

2.328 

Second wing vertical bending 

4 

2. 617 

First fuselage vertical bending 

5 

2. 814 

Outboard pylon lateral bending 

6 

3. 061 

Inboard pylon lateral bending 

7 

' 3. 196 

Third wing vertical bending 

8 

4. 119 

First wing torsion 

9 

5. 013 

Outboard pylon torsion 

10 

5. 155 

Inboard pylon torsion 

11 

5. 387 

Second fuselage vertical bending 

12 

6. 024 

Horizontal stabilizer vertical bending 

13 

6. 7.02 

Second wing chordwise bending 

14 

6. 894 

Fourth wing vertical bending 

15 

7. 891 

Second wing torsion 
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The design model was used as the constant coefficient' example in Section II. 
The state variables chosen were vertical velocity, pitch rate, structural 
mode displacements and rates, control surface displacements, and wind 
model and Kussner function variables. Specifically, the model used in 
Section II consisted of the state vector’ 


^3 5 ^3’ ^ 4 ’ ^4' ^6' ^6’ ^7’ ^7* ^1 1* l'^a’ ^1' ^*2’ ^*3 J ^4' ^*5’ c g^ 

the control vector 


[6 a .6 e ,6 s ] 


the response vector, r, a 14-vector 


r = Hx + Dn 

co n -1 

and the performance index J = E { J r(t) Qr(t)dt}. 

o 

This model was modified to study the effect of parameter variations. To 
reduce the size of the problem the three highest-frequency flexure modes 
(6,7 and 11) were eliminated. Also, the Kussner and wind states (p^, pg, 
p^^P^.Q'g) were eliminated, and the deterministic problem was treated. This 
left an eleven- dimensional state vector 

x T = [x 1 ,x 2 ,...,x n ] T = [w, q,ri r ti 4 , d Q , 6 S ] T 

As initial conditions for this problem we chose x.(O) = y x.. where x.. 

i ii n 

denotes the corresponding diagonal element of the 23rd~order state covari- 
ance matrix corresponding to the optimal control with complete measurement 
used as an initial point for the simplification of Section II. For the parameter 
variations analysis the performance index was chosen to be 

J = f {[x(t)] T Q x(t) + [u(t)] T R u(t)} dt (El) 

o 
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where Q = (H) T QH, R = D T QD and H is obtained from the original H matrix 
by deleting each column corresponding to a deleted state. The^ cross-product 
term 2 x T (H) T QDu was omitted for convenience even though (H) QD was 
nonzero. The parameter variations considered were percentage changes m 
the undamped natural frequencies associated with the three flexure modes. 
Nominal values for these frequencies were assumed to be: 

w '= frequency of first flexure mode = 5.6 rad/sec 

to = frequency of third flexure mode = 14. 5 rad/ sec 

3 

= frequency of fourth flexure mode = 17. 0 rad/ sec 

The vector equation of motion was expressed in the form 

x = Fx + Gu (E2) 


where the following explicit dependence on the u. was assumed: 


43 

-(u/. 

f 65 = » 

f 87 ■ -V 

u 

~ 2 Gi w l' 

f 66 = ” 2 ^3 u 3 

f 88 = " 2 ^4 W 4“ 


PROBLEM DESCRIPTION 

According to theorems 1 and 2 of Section III the optimal insensitive controller 
for a small convex region of admissible system parameters is the optimal 
controller corresponding to a set of parameter values belonging to the bound- 
ary of the convex region which yields the maximum optimal cost. The theory 
gives no quantitative estimate of the possible range of parameters. Similarly, 
the derivations in Appendix D yield qualitative rather than quantitative results 
concerning computational time requirements and approximation accuracy. 
Thus, the major purpose for treating this example was to obtain quantitative 
data from a realistic problem. 
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Two specific questions to be answered with this example were: 

• Is it computationally feasible to compute a second order 
approximation to the optimal performance surface? 

• For how large a range of admissible parameters is such an 
approximation adequate to define an optimal insensitive 
controller? 

Using the adjoint method described in Appendix D, the answer to the first 
question was yes. In answer to the second question the approximation was 
adequate for ± 15%, but inadequate for ± 20% variations from nominal values. 

The optimal cost surface is 

J* = Tr { Px 0 X o T } <E3) 

where P is the matrix function of the parameters which satisfies 

PF 4- F T P - P(GR' i G T )P + Q = 0 (E4) 

The second-order approximation to J* is 

J* = TR [P(u 1 ,to 2 ,u 3 )X o ] + V u J*(Au) + ~ (Aw) T v J*<Aw> (E5) 

where - . -i 
Au^ 

f 

Aw = AWg = the vector of flexure mode frequency variations, 

Aw 4 

w^, Wg, w^ are nominal values of the flexure mode frequencies, 

V J* = ^ - evaluated at o 1( w Q , w,, and v n J* 

dlOg 

is the matrix of second partials of J evaluated at the point w^, w^, w^. 
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The partials of J* were computed from the following equations (which are 
derived in Appendix D): 


2Tr 


i 



r 1 

B J * = Tr J 

3F v 3P + 3F V 3P ’ 

> - Tr J 

!— gr 1gT ll. v ( 

3w.3Wj 

[3co i 3Wj 3w^ SwJ 

I 

IstJj 3 J 


where P satisfies (E-4), V is defined by 

(F-GR -1 G T p)y + V(F-GR _1 G T P) + X Q = O' 


and satisfies 

3 co. 

(f-gr' 1 g t p) t +— (F-GR _1 G T P) + P^- + 

3U| 3w i SUj 

for i = 1, 3, 4. Numerical results from these equations are presented below. 



NUMERICAL RESULTS 


The first and second partials were evaluated to yield the following, approxi- 


mation: 


J* = 21204 + [206, -589, -366] 


Aw. 


Aw. 


Aw. 


+ [Aw^, AtOg, Au^] 

For Aw in the cube Q defined by 

Q = [ Aw : | Aw. | ^ 0. 15 w i , i = 1, 3, 4 }' 


~ -1,45 

-2. 63 

-10. 24 

Au^ 

-2. 63 

10. 66 

15. 46 

^3 

-10. 24 

15.46 

24. 88 

^4 


(E6) 
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the first-order term of (E-6) dominates the second-order term. The point 

/V 

in Q which minimizes J* is readily determined to be 
= 0. 15w^, AUg = -O.lSijg, Aw^ - -0. 15w^ 

The corresponding J* = 24029 compares quite well with the corresponding 
value of J# = 24813. Values of J* and T* at various points in Q are listed 
in Table E2. The accuracy of the estimated value of J* is reasonably good as 
is evident from Table E2. Also J* identifies extremal points of J* correctly. 


Table E2. Optimal Cost and Estimated Cost by Taylor Series 
Expansion at Nominal F 


Parameter Variations 
from Nominal Values (%) 

Optimal 
Cost, J* 

Estimated Cost by 
Taylor Series Expan- 
sion of Nominal F, 

Aw^ 

CO 

3 

<1 

15 

-15 

in 

rH 

1 

24813 

24029 

15 

-15 

15 

21670 

21733 

15 

15 

15 

19608 

19491 

15 

15 

-15 

21089 

21108 

-15 

-15 

15 

21447 

21454 

-15 

-15 

-15 

24091 

23575 

0 

0 

0 

21264 

21264 

-10 

-10 

-10 

22932 

22706 

10 

10 

10 

20088 

19950 

-15 

15 

-15 

20598 

20680 

-15 

15 

15 

19405 

19252 
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To see whether the optimal controller at the maximum cost J^ ax at 15 % Wy 
-15% Un and -15% from the nominal F was indeed insensitive five points in 

O •• f # 

the convex region of the parameter variation were selected, and their optimal 
controllers were obtained. The cost J of the controllers for various points 
in the region were determined. As shown in Table E3, the optimal controller 
at the maximizing point of J* yields lower costs for the points than 
whereas the optimal controller at J* . produces three points whose costs 


exceed 


An attempt was made to perform the same analysis for an enlarged Q varia- 
tion of 20% of nominal values. The resulting optimal cost values and costs 
for the optimal controller at the maximizing points displayed in Table E4 
show that this cube is too large. The optimal controller for the maximizing . 
point yields unstable systems for two vertices of Q. Also, the cost at one 
vertex exceeds the maximum value of J*. This indicates that one would have 
to "back-off" from the maximizing point to derive a satisfactory controller. 
One computationally .feasible way of analyzing this behavior is to analyze the 
locus of roots as the point defining the optimal controller varies. Such an 
analysis was performed for the 15% cube. 


Comparing Table E3 and Figure El, we can conclude that along the diagonal 
line from P5 = (+15% icK, -15% -15% w^) to PI = (-15% w., +15% Ug, +15% 

u .) the optimal cost J* increases. The root loci of the optimal controller 

rx 

of P5 for the points along the diagonal line are given in Figure E2. As seen 
in Figures E2 and E3, if u , w 3 , and u 4 change further than -15% +15% Wg, 

+15% w 4 in the same direction, then the damping ratio decreases and finally 
one actuator pole goes to the right half plane, and hence the system becomes 
unstable. At this point a tradeoff between stability and performance would 
have to be made for the 20% cube. Fortunately, the stability margin need not 
be large since parameter variations have already been included. 
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Table E3. Costs of Controllers at Various Points 



- 



Cost due to the Optimal Controller of 

Parameter Variations 

' 

Awj 15% 

10% (d^ 

5% id 1 

0 % idj 

-15% td 1 

from Nominal Values 

Optimal 
Cost, J* 

Au, -15% u. 

-10% id. 

-5% id. 

0% id. 


15% tdg 
15% u 4 

Aw l 


Ao 4 

A a> 4 -15% 

-10% td 4 

-5% lj 4 

o 

0% td 4 


+15% 

-15% 







• 

-15% 

24813 

24813 

|24888| 

251281 

1255901 


298551 

+ 15% 

-15% 

+15% 

21670 

22455 

22213 

22128 

22146 

22728 

+15% 

+15% 

+15% 

19608 

22036 

20852 

' 20269 

19948 

19641 

+15% 

+15% 

-15% 

21089 

22392 

22200 

22127 

22111 

22577 

+ 10% 
+10% 

+10% 

-10% 

+ 10% 

20088 

21957 

20967 

20500 

20259 

20143 

-10% 

23301 

23375 

23301 

23366 

23570 


256751 | 

+ 5% 

- 5% 

- 5% 

22162 

22478 

22226 

22162 

22219 

23267 

0% 

0% 

0% 

21264 

22072 

21528 

21319 

21265 

21743 

-10% 

-10% 

-10% 

22932 

23118 

22994 

23066 

232 84 

1 25649) 

-15% 

-15% 

+15% 

21447 

22454 

21921 

21767 

21778 

22399 

-15% 

-15% 

-15% 

24091 

24364 

24398 

24655 


25163| 

130799 

-15% 

+15% 

-15% 

20598 

23049 

22309 

22303 

22438 

23268 

-15% 

+ 15% 

+15% 

19405 

24730 

21613 

20374 

19822 

19405 


Note: "] I" means that the value exceeds J* = 24813. 
1 1 max 




Table E4. Costs of the Optimal Controller of (Awi = 20% uj, 
A 103 = -20% ( 03 , Aw 4 = -20% u> 4 ) at Various Points 


Parameter Variations 
from Nominal Values 

Optimal 

Cost due to the Optimal 
Controller of 

( Am ^ — 20 % (i) ^ , AUg — - 20 % (Jg, 

A 1 J 4 = - 20 % 


Au> 3 

Au 4 

Cost, J* 

+ 20 % 

- 20 % 

- 20 % 

26853 

26853 

- 20 % 

- 20 % 

- 20 % 

25566 

| 27110 

+ 20 % 

- 20 % 

+ 20 % 

21756 

23039 

+ 20 % 

+ 20 % 

+ 20 % 

19200 

•25222 

+ 20 % 

+ 20 % 

- 20 % 

21018 

23593 

- 20 % 

- 20 % 

+ 20 % 

21505 

23941 

0 % 

0 % 

0 % 

21204 

23501 

+ 10 % 

+ 10 % 

+ 10 % 

20088 

24346 

+ 5% 

- 5% 

- 5% 

22172 

23146 

- 20 % 

+ 20 % 

- 20 % 

20300 

unstable . 

- 20 % 

+ 20 % 

+ 20 % 

18968 

unstable 


Note: " 


means that the value exceeds ~ 24813. 
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X {+ 15 %w 1/ - 15 %w 3> - 15 %m 4 ) 

O (+10%^ -10 %to 3j -10%o 4 ) 

□ (+57°^ -ST*^ -5 7«w 4 ) 

A C07«w^ 07“ ^ 07 “o 4 5 

V {-157»w^ +157»t03 +157.0)4) 


Figure E2. Root .Loci of the Optimal Controller of 
(Au^ ~ 15 %u)ij AU3 = -15% 103, AW4 - 
-15% <04) for Five Points 



Figure E3. Costs of the Optimal Controller of J^max for 
Five Points Along the Diagonal Line Between 

J *max Point and J *min Point 
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APPENDIX F 

DERIVATION OF THE CHARACTERISTIC POLYNOMIAL 


This appendix derives the characteristic polynomial for an arbitrary m ulti - 
input multi-output feedback control system. We consider plants of the 
following form: 

x = Fx + Gu {F, G) controllable, rank (G) = r 
z = Mx (F, M) observable, rank (M) = m 

where x is an n-dimensional state vector, u an r -dimensional control 
vector and z an m -dimensional vector of measured outputs. Without 
loss of generality, we will assume that the matrix G has a single nonzero 
entry in each column. This will always be the case in aircraft or launch 
booster models where each input drives a separate actuator. Moreover, 
even in cases where G contains additional nonzero entries, we can always 
construct a similarity transformation T such that 

x = Tv 

v = T" 1 F Tv + T" 1 Gu (F2) 

z = HTv 

where the new matrix T ^G has the form 

(F3) 

This is possible whenever rank (G) = r. It can be readily verified by letting 

T = [T^ ; G] , where T^ .is any collection of (n-r) independent columns, 
independent of G. 


T X G = 


FI 



Let the r nonzero entries, of G be denoted by G^ where j-L ?, r, 

and where i takes on r values depending upon j, i. e.j 1 - i(j) - Ul), 
i( 2 ), . . i(r). Then the closed-loop system matrix for a linear fixedt-form 

controller u = Kz will take the following form 


(F + GKM) (l) 


F 


(i) 


F (i) + s G. ;j K^M ( ' e) 
i 


i f i(j) for any j 
i - i(j) 

j=l, 2,... t r 


(F4) 


where the symbol A^ denotes the i row of matrix A. 

In words, the closed loop system matrix is the original matrix P with linear 
combinations of rows of M added to its i(j) ‘ rows.. - 

Now consider the characteristic equation of the closed loop system. This is 
given by 

det (si - F - GKM) = o ( ' F ®‘ 

In order to deal with this equation, we will introduce the following notation 

DET [a (1) , a (2) , ..., a (r) ] = det (A) 


where* 


( (si - F) (l) 

A<*> = 

( -a<*> 


i #i(j) 

i = i(j) for 
3 .= 1 , 2 , . . . , r 


(F6) 


F2 



Further, we will let 0^ = (si - F)^^, 


j — 1, 2j • * • } ^ • 


Verbally, this notation denotes the determinant of the matrix (si - F) 
with rows i(l), i(2), . .., i(r) replaced by -a^, -a^, ... -a 
respectively. Then it follows that 


det (si - F - GKM) = DET [0 (1) + £ G. 


i(l)l K U M 


4<r> + f G i(r)r K r* mU> 1 


U) 


3 * • • 3 


(F 7 ) 


th 


We will now show that equation (F 7 ) corresponds to at most an r order 
polynomial in the gains K. . This will be done iteratively, first for 
r = 1, then r = 2, and finally for general r. 


For r = 1 , note that we are dealing with the determinant of a matrix whose 
th 

i(l) row is a sum of several terms. Using a well known property of 
determinants (see, for example, Hadley, Linear Algebra, page 93 ), this 
gives 

i 

DET [0 (1) + £ G . (1)1 K u M U) ] = DET [0 (1) ] + £ Kj. ^ DET [G. ( 1} 1 M U) ] 

(F8) 


For higher values of r, we simply apply equation (F8) recursively. For 
example, r = 2 yields 
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DET [0 (1) + £ G. (1)1 M U) , 0 (2) + £ G i(2 )2 K 2J 
SL A 


= DET[0 (1) , 0 (2) + SG i(2)2 'K 2i mU) ] + ^ K i i DET[G. (1)1 M U) , 0 ( } 


A 

Ck) 

+ £ G i(2)2 K 2k M ^ 


(F9) 


= DET[0 (1 l 0 (2) J+S K 2x DET[0 (1) , G i(2)2 M (,e) ] 

% 

+ S |DET[G i(1)1 M U) , 0 (2) ] + £ K 2k DET [G. (1)1 M U) , G i(2)2 M (k) ]| 

A 1 

Similarly, for general values of r, we get 
DET [0 (1) + £ G i(1)1 K u M (i) , 0 (r) + £ G i(r ) r K r ^ 


l 


r m 


== DET[0 (1) , '0 (r) ]+ X X 

j=l 1=1 ^ 


.K,.„ DETf0 {I) J . G i<j)j mVA/ - 


U)-, = 

(i) 


<r). 


+ S E X K. K DETfo (1) , ..,G.,. ,, G.,. v • M W2) , 0 (r) ] 

11 22 11 22 

Sum over all nonrepeated pairs of controls 


(F10) 

m m < ) l-r) 

+ 2 • , V; X 1 K * , - K - t DETTm ( 1) G M (%),-^ G i(i ) 3 M^.00 ] 

V 1 3 i*i 3 p a p DETC0 ’••' G i(3 1 >;5 1 M Vp 

Sum over all nonrepeated .groups of p controls (rl / p! (r-p)I terms) 


m m 


+ E . E . . K u _ ...K rA DE.TtG 1(1)1 G 1(r)r 


li = 1 « r =1 


Ui> n M Ur) 1 
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Interpretations 


The various terms of the polynomial (F 10 ) have the following interpretations 

1) DET^f^, . . . , 0^} is the n th order characteristic 
polynomial, D(s), of the open loop system. 


2 ) DET(0 (1) G. (j) . M U> , . . ., 0 (r) ) with G i( . } . M U) in the j th 

position is the (n-1) order (or lower) numerator polynomial 
of the (j th input) -to -{£ th output) transfer function. This can be 
verified by computing the transfer function: 

0 


— (s) = 
u. 

3 


-M U) (sI - F)" 1 j G.,. v . 

! 1(3)3 ; 

» 1 


n 


1 th. 

= ^(i) k ^~ m ik Ci(j) k cofactor of si - F] G.^.. 


thr DET [<4<1> mU> « s<r);] G i<j)j 


.... G i(j)j M U) , .... 0 <r) ] 


The numerator polynomials will be denoted by N (s). 

Jv 
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3 ) 


For higher-order terms of equation (F10) the DET [ ] expressions 
are so-called "coupling numerators" [18], which appear whenever 
several feedback loops exist simultaneously. We will denote these 


hh' 


J P 


*1*2' * **p 


polynomials by N (s), i. e., 

N* 1 P (s) = DET [0 (1) , G.,, „ M V *P', . . , 0 Vi/ ] 


h m * *^> 


i(lp)lp 


rUr.) »■ 

* 

(FI 1) 


Simplifications 


Equation (F10) can be significantly simplified by using the following two properties. 


hh 9 * ‘ 3 p 

1) N (s) = 0 

^1^2* * * ^p 


whenever two or more elements, of the sequence Ay -^p 

are identical. This is because the determinant of a matrix with 
two proportional rows or columns is zero. 


2 ) 


hh' * ‘^p hh' * 

N (s) = (+) N ^ 

' e i i 2‘ ' * x p *1*2* * 



whenever the sequence A A, . .., A' is an arbitrary permutation 

1 6 P 

of the sequence A*, A 0 , A . The correct sign is obtained by 

noting whether the permutation is even (plus sign) or odd (minus 
sign). This follows from the fact that an interchange of any two 
rows or columns of a matrix only changes the sign of its determinant. 
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Using property (1), we can readily conclude that the maximum order in K of 
the polynomial (P10) is not equal to the number of controls (r) but equals 
min (r, m)-. To show this,- just compute a higher order term, i. e. , let 
r p > m. Then definition (Fll) for 




requires that .we substitute some rows of M into (si - F) 


more than once, so the polynomial must be zero. We can thus eliminate 
a whole host of terms in equation (F10). 


Still more terms can be eliminated by using property (2). Note, in particular, 
that the p th order summations 


m m 

E • ♦ • E 

|,=1 I =1 

*1 *p 

range over all possible combinations of p numbers out of m, including groups 
with repetitions of numbers and groups which are permutations of other groups. 
The former will yield vanishing polynomials while the latter's polynomials 
will differ only in sign from other polynomials. We can thus restrict the 
summation to nonrepeated groups only, and we can collect all groups which 
are permutations of each other to form a single multiplier of the polynomial 
corresponding to, say, the naturally ordered permutation. This gives the 
following reduced characteristic equation: 


min(r, m)\ 
P(s) = D(s)+ E 
P=1 


h m * ' 2 p 


E E N 

3 r ; * 3 p V h 

♦ * 


. (s.)[e(±)K ...K ] 

...£ a 3 11 Vp 

p t 


(F12) 


sum over all permutations 
of , 

sign 


of A 1 . . . X with appropriate 


sum over all naturally ordered nonrepeated 
groups of p measurements out of m 


sum over all naturally ordered 
nonrepeated groups of p controls 
out of r 
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APPENDIX G 

POLE-PLACEMENT EXAMPLES 


This appendix describes two pole-placement problems which were used to 
debug, exercise, and demonstrate the pole-placement algorithm developed 
in Section IV. The first problem involves rigid-body controller designs for 
a sixth-order lateral- axis model of the F4 aircraft (Ref. 6), with rigid -body 
pole placement taken as the design objective. The second problem is con- 
cerned with flexure control of a 17th-order pitch-axis model of the C-5A 
transport. In this case, the design objective is taken to be placement of 
rigid body and certain critical flexure poles. These examples serve two 
functions — (1) they show that the algorithm works, and (2) they illustrate 
how it can be utilized in control design and sensor choice problems. 


THE F4 LATERAL-AXIS PROBLEM 

Disregarding servos, flexure, and other high-frequency phenomena, the 
lateral axes of the F4 can be modeled by the following equations: 

» 

x = Fx + Gu 

where 

__p 

r 

0 

6 r 
6 as 


- roll rate (stability axes) 

- yaw rate (stability axes) 

- angle of sideslip 

- angle of bank 

- rudder actuator deflection 

- aileron/ spoiler actuator deflection 
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u = 


rudder actuator input 


u. 


R 


u 


AS 


aileron/ spoiler actuator input 


and where the matrices F and G take the following values for a flight condi- 
tion at Mach = 0. 5, altitude = 5000 feet: 


-1.7680 

0.4125 

-14. 520 

0. 0000 

-0. 0007 

-0. 3831 

6. 038 

0. 0000 

0. 0016 

-0. 9975 

-0. 155 

0. 0586 

1. 0000 

0. 0000 

0. 000 

0. 0000 

0. 0000 

0. 0000 

0. 000 

0. 0000 

0.0000 

0. 0000 

0. 000 

0.0000 



* 


“ 


2. 031 

8. 9520 


0 

0 

-3. 398 

-0. 3075 


0 

0 

0. 028 

-0. 0036 

G = 

0 

0 

0. 000 

0. 0000 


0 

0 

20. 000 

0. 0000 


1 

0 

0. 000 

-10. 0000 


0 

1 


To represent unity gain actuators, the nonzero values in G (i. e. , Gg^ and 
G ) are equal to 20 and 10, respectively. For our own convenience, 
however, we will consider these values to be unity and scale down all final 
gains by a factor of 20 or 10, as required. 


The dynamics of the above system are characterized by three principal 
modes; 

(1) The spiral mode, corresponding to a root close to the origin 
(-0. 0156) 

(2) The roll subsidence mode, corresponding to a root at -1. 85 

(3) The dutch-roll mode, corresponding to a complex conjugate 
pair of roots at - 0. 219 ± j 2. 48 

In addition, there are two actuator modes, corresponding to roots at -20 
and -10, respectively. 
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The "quality" of the lateral-axis dynamics is judged to a large degree by the 
root locations associated with the three principal modes. In particular, the 
spiral root should correspond closely to a pure integration; as it does in the' 
free aircraft above. This allows pilots to hold nonzero bank angles and turn 
rates without stick commands. 

The roll subsidence root should be fairly large, since it defines the rate of 
response of roll rate due to aileron inputs. Values of -3 or -4, for example, 
are much preferred over the rather sluggish -1. 85 value. Finally, the dutch- 
roll roots should be damped at least by a ratio of 0. 25 and should exhibit 
natural frequencies in the vicinity of 2. 5 rad/sec. This assures tolerable 
sideslip and/or roll-rate oscillations in response to lateral control inputs or 
disturbances. Note that the free airframe has good dutch-roll frequency 
characteristics but falls short of the damping requirement. 

Given these requirements on pole locations, the lateral axis problem makes 
an attractive pole placement example. In the paragraphs below, we consider 
the possibility of satisfying the requirements with several different sensor 
complements. In each case, the algorithm of Section IV was used to assess 
pole placement capability and to compute controller gains. 


Case 1. Four Measurements 


The following measurements are usually available on the aircraft: 


y = Mx = 


se 


se 


0 


- sensed roll rate (body axis roll rate) 

- sensed yaw rate (body axis yaw rate) 

- lateral acceleration at the accelerometer station 

- bank angle 
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The measurement matrix corresponding to these signals is given below. 


0. 9985 

-0. 0541 

0 . 00 

0 

0 . 00 

0 . 000 

0. 0541 

0. 9985 

0 . 00 

0 

0 . 00 

0 . 000 

0. 6084 

-2. 3640 

-27. 33 

0 

-17. 85 

-3. 695 

0.0000 

0 . 0000 

' 0 . 00 

1 

0 . 00 

0 . 000 
J 


The various terms of the closed-loop characteristic polynomial [equation (44), 
main text] for this set of measurements were computed with the algorithm 
and are summarized in Table Gl. This table consists of 15 nonzero coeffi- 
cient vectors — one zero-order term, D, eight first-order terms of the 
form and six second-order terms of the form Each vector is 

identified in Table Gl according to its order and the control sequence 

, 3 ' } and measurement sequence . . . , use(3 to generate it. 
From these sequences, the gain multiplier, ^(K), corresponding to each 
vector can be reconstructed as follows: 


General multiplier = £ (±) K. . K. . ... . K- . 

3 1*1 3 2*2 3 p*p 

summed over all permutations of 

• • • > *p^ 

Specific multiplier (p=l) = K. . 

3 rr 

Specific multiplier (p=2) = (K. . K. . - K. K. ) 

3 jl *! ^ 2*2 3 1*2 3 2*1 


In a neighborhood of the initial gain K Q = 0, -the number of arbitrary poles 
p is given by the rank of .a matrix formed from all the first order terms 
of Table Gl. This rank turns out to be six, so that ,all of the system poles 
can be placed' arbitrarily. Consequently, the .following closed-loop pole 
assignments were made.: 
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Table Gl. Polynomial Terms for F4, Four Measurements 


Control Sequence . . j K ] 

Measurement Sequence 
Coefficient Vector (n,n 0# . .n_)T 

i Ci Q 


0 



• 3 fcoo£ 

02 

•2329E 04 

•176oE 04 

•6846E o3 

•2762E 03 

•3231E o2 

1 

1 

1 

-•1175E 

01 

.3706E 03 

•408 S E 02 

-.2174E o2 

••2212E 01 

•OOOOE 00 

1 

1 

2 

•81-S9E 

02 

•2915E 02 

• 6665 E 02 

•3923E oa 

•3283E 01 

•OOOOE 00 

1 

1 

3 

•»29S1E 

02 

•3957E 04 

•2454E 04 

•532SE o3 

•2112E 03 

•1785E 02 

1 

1 

4 

.37 l7 E 

03 

•4436E 02 

»»195gE 02 

*•2031 E 01 

. 0000 E 00 

• 0000 E 00 

1 

2 

1 

•3U3E 

01 

••9982E 03 

••1453E 03 

-.1S39E 03 

••8956E 01 

•OOOOE 00 

1 

2 

2 

-•58532 

02 

••5607E 02 

«458qE. ol 

-.3179E 5l 

-•l770E 00 

•OOOOE 00 

l 

2 

3 

•25&5E 

03 

•8493E 03 

•647fcE 03 

•7529E 02 

. 7615 E 02 

• 3695E oi 

l 

2 

% 

••9997E 

03 

-.1449E 03 

• • 183sE 03 

-•8952E 6l 

•OOOOE 00 

•OO 00 E 00 

2 

1 2 
1 2 

• 4353 E -10 

• 3351 E 01 

•2979E 02 

•OOOOE O0 

•0000 E 00 

• 0000 E 00 

2 

1 2 
1 3 

-•58J9E 

01 

•1834E 04 

•2433E 02 

•1517E o3 

•OOOOE 00 

•OOOOE 00 

2 

1 2 
1 4 

•18j2E 

00 

•1611E 01 

•OOOOE 00 

• 0000 E 00 

• 0000 E 00 

•OOOOE 00 

2 

1 2 
2 3 

•1075E 

03 

•9931E 02 

••4928E ol 

• 1529E ,o2 

•OOOOE 00 

•OOOOE 00 

2 

1 2 
2 4 

-•3346E 

01 

-•2975E 02 

•0000E 00 

•OOOOE 00 

•OOOOE 00 

•OOOOE 00 

2 

1 2 
3 4 

-•I836E 

04- 

* *24032 02 

- • 1523E 03 

•0000E 00 

•OOOOE 00 

♦OOOOE 00 
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Spiral (s + 0) 

Roll subsidence (s + 4. 0) 

Dutch roll (s 2 ' +1.25 s + 6.25) 

o 

Actuator (s + 30 s + 450) 

and the algorithms executed seven Newton-Raphson iterations before the con- 
vergence test was satisfied. The final gains were 

-41.03 110.60 0.1649 -28.64 

17.24 -46.65 0.0000 0.00 

The convergence test required that the spiral root be located to an accuracy 
of ± 0. 001, the roll subsidence root to an accuracy of ± 0. 01, dutch- roll coef- 
icients to ± 0. 01 and the actuator coefficients to ± 1.0. An independent check 
of the closed-loop roots showed that all conditions were satisfied. 

It is important to recognize that each iteration of the algorithm requires a 
square [in this case (6x6)] Jacobian matrix. Such a matrix was obtained from 
our original (6x8) Jacobian by choosing the first n independent columns and 
deleting the rest. This arbitrarily removes two degrees of freedom and 
accounts for the fact that the final values of K 23 and K 24 are zero. 

Case 2. Three Measurements 

ti 

Having shown that four measurements are sufficient for arbitrary pole place- 
ment, we next proceeded to delete bank angle from the measurement comple- 
ment. The program was again used' to compute coefficient vectors, to assess 
rank, and to carry out Newton-Raphson iterations. 
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A bit of reflection will verify that the set of coefficients vector for three 
measurements must consist of all vectors in Table G1 which do not involve 
the fourth measurement. This gives ten vectors -- one zero-order term, 
six first-order terms and three second-order ones. The rank of the first- 
order vectors is again six, so all system poles can again be assigned arbi- 
trarily. For the same assignments used above, the algorithm proceeded to 
execute six iterations before terminating with the following gains: 


K 


2.39 -10.44 -1.19 

-16.67 83.57 6.56 


The desired root accuracies were again verified by an independent check of 
the closed-loop roots. 


Case 3. Two Measurements 


Knowing that three measurements is enough, one is tempted to try two -- 
so we deleted the lateral acceleration signal as well as bank angle. This 
leaves only two rate gyros as the sensor complement. 

The set of coefficient vectors is again a subset of Table Gl, found by re- 
moving all vectors involving measurements 3 or 4. This leaves six vectors 
-- one zero-order, four first-order, and one second-order. The rank of the 
first-order terms is now four, so that only partial pole placement is possible. 
The initial gain K q = 0 was again selected and the spiral, roll subsidence, anc 
dutch- roll roots were defined as above. This exhausts the four roots which 
may be arbitrarily specified. We simply have to accept whatever the algor- 
ithm gives us for the remaining roots. Final gains and unspecified roots 
were found in three iterations: 



6(s) = fi 2 ■+ 27. 06 s + 149. 1 

actuator roots at -19. 35, -7.71 

We see that even though the actuator poles could not be' placed arbitrarily, 
they take on quite reasonable values. 

It appears, therefore, that only two measurements can suffice to achieve 
modal control of the lateral axes. 


THE C-5A FLEXURE PROBLEM 

Appendix iE contains descriptions of a 23rd-order flexure model for the C-5A 
pitch axis. This model was utilized here for high- order pole placement 
studies. However, since gust dynamics and Kussner lift growth states do not 
enter the feedback loop, these six degrees of freedom were deleted. This 
leaves a 17th-order model, comprised of two rigid body states, two states 
each for the 1st, 3rd, 4th, .6th, 7th and 11th flexure modes, and one state 
each for symmetric aileron, elevator, and spoiler actuators. 

Open- loop .poles for the 17th-order model are given .in Table G2 together 
with a set of ".desirable locations" abstracted from quadratic optimal flexure 
control designs and from the LAMS controller (Ref. 10). These locations do 
not correspond to any one controller but were chosen to reflect general 
trends exhibited by several designs — namely, that flexure controllers move 
rigid-body poles, increase damping of modes 1, 4, and 7, and leave the 
remaining flexure poles pretty well alone. This same trend was taken to be 
the design objective for the pole- placement computations reported here. 
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Table G2. Open-Loop and Desirable C-5A Roots 


Mode 

Open-Loop Roots 

Desirable Roots 

Rigid Body 
1st Flexure 
3rd Flexure 
4th Flexure 
6th Flexure 
7th Flexure 
11th Flexure 

-1.03 ±j 1.31 
-0/901 ±j 5.57 ,! 

-0.462 ± j 14. 5 
-1.05 ±j 17.0 
-0. 430 ± j 19. 3 
-0. 960 ± j 20. 3 
-4. 29 ± j 35. 7 

-1.07 ±j 1.80 

-4. 00 ±j 6,00 

Unchanged 

-4. 50 ± j 18.0 

Unchanged 

-11.0±j 22.0 

Unchanged 


A sensor complement composed of two accelerometers, one rate gyro, and 
one surface position transducer was used in the computations. The acceler- 
ometers were located on the wing, one near midspan and the other near the 
wingtip, while the rate gyro was located near the tail. These locations corre- 
spond to the LAMS controller (Ref. 10). The position transducer measured 
aileron actuator deflection. The measured sitnals are therefore summarized 

- midspan acceleration 

- wing tip acceleration 

- pitch rate at tail station* 

- aileron actuator deflection 

Each signal is a linear combination .of rigid body states, bending mode states, 
and actuator states. 



■ *Due to an error in programming this signal consisted of rigid-body 
pitch rate plus mode slopes times mode positions rather than mode 
rates times mode slopes. 
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In a neighborhood of the open- loop system, these measurements yield 


max 


= rank 


K q =0 


N 


N 4 N 1 • ‘ 


n: 


~ 8 


This is indicated as an approximate value because the numerical determina 
tion of the rank of high order matrices is subject-to uncertainty. . In the pres- 
ent case, the matrix of numerator coefficients spanned 8-dimensional space 
with certainty but possibly also 9-dimensional space. Components in the 9th 
coordinate direction, however, were near the least-bit level of the machine 
and thus of little practical utility. At any rate, eight closed-loop poles could 
be assigned arbitrarily in a neighborhood of the open- loop system. 


Given this - information,, three separate runs were made with the initial condi- 
tion K = 0, The first run placed four poles (rigid body plus 1st flexure), the 

o i . , 

second six -poles (rigid body plus 1st and 4th flexure), and the third placed 
eight poles (rigid body plus 1st, 4th and 7th flexure). Results are summarized 
in'Table G3. In each case, open- loop actuators at -100 (elevator),’ -10 
(spoiler), and -6 (aileron) were assumed, and a square Jacobian matrix was 
obtained by allowing only selected feedback gains to vary. The resulting gain 
structure as well as the number of Newton-Raphson iterations are also shown 
in Table G3; 


As the table indicates, each run achieved its objective of moving a specified 
group of poles into the desirable locations given in Table G2. However, some 
of the remaining unassigned poles were destabilized in the process. In Run 1 
the damping of 3rd, 4th and 7th flexure modes was. reduced; in Run 2 the 7th 
mode was actually driven unstable; and in Run 3 an actuator root was driven 
unstable. Fortunately, none, of these runs fully exhaust the pole placement 
capacity offered by the sensors. Away from the origin, the previously weak 
9th order of' rank strengthened sufficiently >to yield 


max 


-= rank": 


EN. 


K o K run 3 


5K T 


= 9, 
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Table G3. Pole Placement Runs for C-5A 


Mode 

Rigid Body 
1st Flexure 
3rd Flexure 
4th Flexure 
6 th Flexure 
- 7th Flexure 
11th Flexure 


Run 1 

4 Poles Placed 

-1. 07 ± j 1.80 

-3. 95 ± j 6.00 

-0. 253 ±j 8.68 

-0. 182 ±j 16.5 
-0. 258 ± j 19.3 
-1. 59 ±j 18.6 


-4. 27 ± j 35.7 


Run 2 

6 Poles Placed 


-1. 07 ± j 

1 . 80 

-4. 00 ± 3 

6 . 03 

-1. 07 ± j 

14. 9 

-4.49 ± 3 

18.3 


-0.448 ± j 19.4 
+0. 004 ± j 17, 8 
-4. 54 ±j 36.8 



Gain 

Structure 


Elevator 

Spoiler 


Z 1 

Z 2 

z 3 

Z 4 

X 

X 

X 

X 




x x 

X X 
X X 


Run 3 

8 Poles Placed 

-1.07 ±j 1.80 

-4. 09±j 6,07 

-0. 955 ±j 16.4 
-4. 59 ± j 18.0 
-0.357 ±j 19.3 
tIO, 9 ± j , 22. 2 
-4.08„±j 35.0 


-103. 0 
+ 6 . 15 

- 10. 1 


X 

X 


X 

X 

X 

X 


X 

X 




Run 4 

9 Poles Placed 

-1. 06 ± j 1 . 79 
-4. 00 ± j 5. 94 
-0. 506 ± j 15. 8 
-4. 53 ±j 18.3 
-0. 265 ±j 19.3 
-10. 9 ± j 22, 0 
-4. 15 ±3 35.1 


-108. 0 

- 0. 972 

- 2. 52 


x x 

XXX 
X X 


Total number of iterations in 20-step incremental stepping procedure 












which means that one additional pole could be placed arbitrarily. This extra 
degree of freedom was used in Run 4 to move the unstable actuator pole back 
into the left-hand plane. The resulting final pole constellation (Table G3) 
shows that our design objectives can be satisfied by utilizing the full capacity 
of the sensor complement. All desired locations are achieved and the re- 
maining poles are approximately unchanged. Gain magnitudes associated 
with Run 4 compare favorably with the LAMS controller. 


An interesting byproduct of these pole-placement computations is what might 
be called a generalized root locus plot. Namely, if an incremental stepping 
procedure is used to move from the open-loop poles to the desired closed- 

loop poles, i. e. 


Ms, a) “ <2 

a = 0, Act, 2 Act, 


\ . ,(s) - k (s) 

desired open 


loop 


+ ^ (s) 

open 

loop 


then pole positions as a function of or form a familiar single-parameter root 
locus which shows how unspecified roots migrate as the assigned roots move 
toward their specified positions. Such plots are shown in Figures G1 and G2 
for Runs 1 and 2 of Table G3. Both figures indicate that acceptable levels of 
damping can be obtained for the unspecified roots if a value of o' somewhat 
less than unity is accepted for the assigned roots. 
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Figure Gl. Root Locus: Run 1 
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*?Eigur.e.-'G2,. .Rpot *Locus: 'Run '2 
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APPENDIX H 

PROCEDURES FOR COMPUTING COEFFICIENT VECTORS 


This appendix describes two methods for computing coefficients of polynomials 
appearing in the general characteristic equation of Appendix F. These parti- 
cular computations comprise about one half of the entire computational load 
associated with the pole placement algorithm of. Section IV. 


METHOD 1 

This procedure is presently used in the algorithm. It mechanizes the computa- 
tional method of Reference 19 for solving the general eigenproblem det (sB+A), 
where A and B are arbitrary matrices, possibly singular. (The DET [. . . ] 
expressions of Appendix F are in this form. ) The method consists of reducing 
the original problem to a lower -order nonsingular one by repeated Gauss 
reductions of matrices B and A. The nonsingular problem can then be solved 
by standard methods. The reduction goes like this: 


• Via Gaussian elimination, reduce B to 

l 


B' = 


B ll!’ B 12 


0 


0 


where is upper triangular. 

AH elementary operations performed on B are also performed on A, which 
gives 


A' = 


A ii! A i 2 ’ 

A 2l]" A 22~_ 


HI 



Via Gaussian elimination, reduce' A' to 


A" = 


"A" I' A" 
A n i A i2 


. O' 


A" 

■I A 22. 


where A'' is. upper - triangular; All' elementary operations performed on 

j 

A 1 are; also performed, on B' 1 ,. giving 


B" = 


B ir j 

< 

1 

0 


We now have- 


det .(sBTA), =■ ±. det < s' 



B 12 _ 

• 

A" 

,11 

A" 

12 



+ 



■ 0 - 

'o 


,0 

A-"' 1 
22_ 


=-±. det A” 2 . detCsB'.^ + A'^ 1 > ; 


where the. sign is. determined 1 by r.oW and column operations, performed 
in. the two- Gaussian elimination steps.. 


If is nonsingular;, the problem! has: been, reduced 1 to a- standard one-. 
Otherwise,, we- let, B. - B^,, A =• A'^. and return to Step (1). 


The remaining nonsihgular ei'genproblem is presently being solved by computing 
eigenvalues: via. QR'. transformations: and: reconstructing- coefficients from these. 
Total SDS 9300- computings times, are- 3 - to h seconds, per 17th-order’ coefficient 
vector., 
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METHOD 2 


This method utilizes the fact that higher-order generalized n um erators can 
be expressed in terms of first-order ones. This offers significant com- 
putational savings when repeated computations must be performed for • 
different measurement matrices, M(fi, y). 


Expressions for Higher Order Numerators 


An alternate way to derive the characteristic equation of Appendix F is via 
the transfer function matrix -H(s) relating measurements z and control 
inputs u; i. e. , 

z{s) = [-H(s)Ju(s) 

u(s) = Kz(s) + u (s) 
c 

[1 + H(s)K]z(s) = [-H(s)]u (s) 

P(s) = D(s) det (1 + H(s)K) 

Expanding this polynomial and comparing the result term for term with 
equation (F12) of Appendix F provides expressions for the higher order 
numerators. This process is carried out for a 3 -output, r-input system 
below, from which the general relationships can be deduced. 
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P(s) = D(s) det 


i + SH llKll 

EHi.K i2 

EHiiK i3 


SH2i K il 

1 + EH 2i K 12 

SH 2i K i3 . 


EH S1 K 11 

EH gi K i2 

1 + eh 31 k. 3 


S u + S 22 + S 33 + 

S 11 S 22 ' S 12 S 21 + E 11 S 33 

S 13 S 31 + ^22^33 


+ S ll 2 22 2] 33 " 2 11 S 23 S 32 “ 2 12 E 21 2 33 + S 12 S 23 S 31 
+ S 13 S 21 S 32 ~ I ’l3 S 22 S 31-' 

where E^ = S 

Note that the terras which are first-order in £ may be written as 

v 

S 11 + S 22 + S 33 ~ ^ =1 ? =1 Hkftk.” * \ 

K h i 

Similarly the second-order terms may be expressed as: 

y. WWW* 

J l 3 2 

H . Ho- (K. ,K o - K gK. J + 

h 1 3 3 2 h 1 h 3 . h V 


H 9 . Hn- (K. gK. o “ K. oK- 2) 
2 h 33 2 3l 2 ^2 3 h 
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This expression can be simplified by observing that terms are skew- 
symmetric in j 2 , i. e. , interchanging the roles of and j 2 in a term 
merely changes the sign of the term and in particular if 3 ^ = j 2 then the 
term is zero. Thus the second-order terms are given by: 


2 2 (H H ~ H H )(K .K 2 - K „K .) + 

j 2 >j x lj i 2j 2 x 3 2 2 h h 1 h 2 h 2 h 1 


(EL . H„. - H . H q . )(K . .K. o~K. r,K. .) + 

X 3l 3 3 2 J 3 2 3 3 i h 1 3 2 3 3i 3 3 2 X 




H„, H„. )(K 9 K. o-K. ,K . J 
2 h 3 3l 3 i 2 3 2 3 3 X 3 J 2 2 


or, in general form 


2 _ 2 
3i ^1 


CE (±) H H ] [£ <±> K K ] 

* 1^1 *2 j 2 P j 1*1 j 2*2 


I L 


1: 


1 J 2 *12 

a" T 

all permutations of 

permutations of 

^—all ordered nonrepeated groups of 2 measurements out of 3 

i 

L all ordered nonrepeated groups of 2 controls out of r 


Similar expansion and reduction of the third -order terms yields 


4,3 Ws ^ (±> VW % (±) Vx'WW 


From these expressions, it follows that the general closed-loop polynomial 
for r inputs and m outputs is given by 
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( min(r, m) 

P(s) = D(s> 1+ S . S , . ~ p ..j, 

( k-1 3^*-! * ^1* * ' ^ ^ ^ 


E [E (±)H i,j," H i, i, 1 


I 


[S (±)K. . ..K . ] 

P ^1*1 ^k^k ) 

0 . ' 


and comparing with equation (FI 2) of Appendix F yields 


n/ .Ns) = D(s) [S (±)H A1 . 

P. *L 3 1 2 ] 2 

3 


H> 


"k J k 


] 


1 


D(s) 


k-1 


j 1 J2 ^k 
[E (±) N. Ni ..N 3 
p. *, * \ 


This expression defines all higher -order numerators in terms of first-order 
ones. Note, however, that the definition involves a sum of (nk) order 
polynomials divided by the n(k-l) th order polynomial D k We know from 
Appendix F that each higher -order numerator is, in fact, an n -order 
polynomial (or lower). This means that D must be a perfect factor of 
the sum of (nk) th order polynomials - a fact which could lead to numerical 
problems in computing N^U* * * ^ from the first-order terms.. The following 

1* * * K 

procedure is suggested: 

• Multiply the definition by D k_1 and express the result in 
coefficient -vector -form: 


- ^1* * * ^k 



N 
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where D is an (nk x n) -dimensional matrix constructed from the 
coefficients of D(s) and N is the coefficient vector of 


h 


2 (±) N . . . N 


P. 

3 


j k 


• Solve for N. . via pseudo inversion, i. e. , 

~ i 


3i • • - J%. 

1 ' K -T _ _1 _rp 

N, , = (D 1 D) i D 1 N 

1* * * lc 


This should minimize numerical problems associated with the polynomial 
division. Note that the quantity (D D) 3D needs to be computed only once. 
It can be stored and reused for computation of all k^-order numerators 
xor any given set of first-order numerators. 


Computation of First-Order Numerators 


All that remains now is to compute the first-order numerators. This can 
be done via Leverrier's algorithm (also called Fadeeva's methodi) which 
computes the determinant and all of the cofactors of (sI-F)" 1 . That is. 


tZadeh, L.A. , andC.A. Desoer, Linear System Theory, McGraw-Hill, 
New York, 1963. 



(sI-F) 


-1 


Bs nl +B - s a 2 + . . . + B, 
B(s) n n-1 1 


D < s > s n + a s 11 + d ,s n - 2 +...+d, 
n n-1 1 


B = I 
n 


d = -Trace (B F) 
n n 


B . = B F + d I 

n-1 n n 


d , = -l/2Trace (B , F) 

n-1 ' n-1 


B k “ B k+1 F + ^l 1 


n+l-k 


Trace (B k F) 


B 1 = B 2 F + d g I 


d^ TracetB.^F) 

n 


0 = BjF + d x I 


Once all the cofactors are known, individual first-order, numerators are 
obtained from the relation 


N i< s > - * - m & B k«3) <s > G i(:)j 

This expression shows that repeated runs with different M matrices are very 

-T - -1 -T 

economical, since the B(s) matrix {and also the (D D) D expressions above) 
need not be recomputed for a new M matrix. 
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